ISSN 0022-5096 


/ YORK 


¥ 


* TORONTO 





JOURNAL OF THE MECHANICS AND PHYSICS OF SOLIDS 


Editor-in-Chief 
J. R. WILLIS 


School of Mathematical Sciences, University of Bath, 
Claverton Down, Bath BA2 7AY 


Editorial Advisers 


R. HILL (Editor-in-Chief 1952-68)—-Cambridge 


P. CHADWICK—East Anglia R. J. CLIFTON—Providence 
K. L. JOHNSON—Cambridge J.W. HUTCHINSON—Cambridge (MA) 
J. R. RICE—Cambridge (MA) 


Publishing offices: Headington Hill Hall, Oxford OX3 OBW, U.K. (Oxford 64881) 


Annual Subscription Rates. 1990 (including postage and insurance ) 


Annual institutional subscription rate (1990) DM 660.00. 2 year institutional rate (1990/91) DM 1254.00. 
Personal subscription rate for those whose library subscribes at the regular rate (1990) DM 189.00 


Issued six issues per annum. One volume per annum. 


Subscription enquiries from customers in North America should be sent to: Pergamon Press Inc., 
Maxwell House, Fairview Park, Elmsford, NY 10523, U.S.A., and for the remainder of the world 
0: Pergamon Press plc, Headington Hill Hall, Oxfor dOX3 0BW, U.K. 


Microform Subscriptions and Back Issues 


Back issues of all previously published volumes, in both hard copy and on microform, are available direct from 
Pergamon Press offices. 


Copyright © 1989 Pergamon Press plc 


It is a condition of publication that manuscripts submitted to this journal have not been published and will not be simultaneously submitted 

or published elsewhere. By submitting a manuscript, the authors agree that the copyright for their article is transferred to the publisher if and 

when the article is accepted for publication. However, assignment of copyright is not required from authors who work for organizations 

which do not permit such assignment. The copyright covers the exclusive rights to reproduce and distribute the article, including reprints, 

photographic reproductions, microform or any other reproductions of similar nature and translations. No part of this publication may be 

reproduced, stored in a retrieval system or transmitted in any form or by any means, electronic, electrostatic, magnetic tape, mechanical, 
photocopying, recording or otherwise, without permission in writing from the copyright holder 


U.S. Copyright Law Applicable Users in the U.S.A. 


Photocopying information for users in the U.S.A. The Item-Fee Code for this publication indicates that authorization to photocopy items for 
internal or personal use is granted by the copyright holder for libraries and other users registered with the Copyright Clearance Center (CCC) 
Transactional Reporting Service provided the stated fee for copying, beyond that permitted by Section 107 or 108 of the United States 
Copyright Law, is paid. The appropriate remittance of $3.00 per copy per article is paid directly to the Copyright Clearance Center Inc., 27 
Congress Street, Salem, MA 01970 

Permission for other use. The copyright owner's consent does not extend to copying for general distribution, for promotion, for creating new 
works or for resale. Specific written permission must be obtained from the publisher for such copying 

eo” The text paper used in this publication meets the minimum requirements of the American National Standard for Information Sciences 
Permanence of Paper for Printed Library Materials, ANSI Z39.48 1984 


The Item-Fee Code for this publication is: 0022. 5096/90 $3.00 + 0.00 








J. Mech. Phys. Solids Vol. 38, No. 1, pp. 1-25, 1990. 0022-5096/90 $3.00 + 0.00 
Printed in Great Britain. © 1989 Pergamon Press plc 


PLANE STRAIN CHARACTERISTICS THEORY FOR 
SOILS AND GRANULAR MATERIALS WITH DENSITY 
DEPENDENT YIELD CRITERIA 


I. F. COLLINS 


Department of Theoretical and Applied Mechanics, School of Engineering, University of Auckland, 
Auckland, New Zealand 


(Received 2 August 1988 ; in revised form 6 February 1989) 


ABSTRACT 


AN ANALYSIS of the quasi-static, plane strain, deformations of a rigid/plastic material whose yield stress 
depends on the current material density is presented. Both steady and unsteady flows are studied. Such 
constitutive relations are widely used to model the yielding behaviour of soils and granular media in what 
is generally known as “‘critical state soil mechanics’. It is shown that, contrary to a widely held belief, the 
stress and velocity characteristics for such materials do coincide. The differential relations along the 
characteristics are derived and shown to have a structure typical of hardening materials. In addition it is 
shown that a second family of weak discontinuities, embedded in the material, can exist and that they give 
rise to the familiar ‘“‘“Sokolovskii stress characteristics” in the critical state limit. 


INTRODUCTION 


As IS WELL known, the classical theory of slipline fields developed to describe the 
plane strain deformation of metals does not generalize readily to granular materials. 
The simplest models of the deformation of soils and granular media are based on the 
Mohr-—Coulomb failure criterion, which states that the material fails when the shear 
stress t acting on the surface of a material element attains the critical value: 


lt] = c+otany, (1.1) 


where o is the compressive normal stress and c the cohesion and wy the angle of 
internal friction are material constants. When expressed in terms of the in-plane stress 
invariants, this condition becomes 


g = psinw+ccosy, (1.2) 


where 


p=3(0;;+022) and gq =[4(01;—022)? +072] (1.3) 


are the mean pressure and maximum shear stress respectively. We are employing 
the usual convention in geomechanics of regarding compressive normal stresses as 
positive. 
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This failure condition, together with the two quasi-static equilibrium equations, 
forms a statically admissible system of three equations for the three in-plane stress 
components. This system is quasilinear and hyperbolic with the two families of 
characteristic curves inclined at angles of +(2/4—w/2) to the major (more com- 
pressive) principal stress direction. The Coulomb failure condition (1.1) is actually 
attained by the traction components on these characteristic curves. A comprehensive 
account of this theory is give by SOKOLOvsKiI (1965) together with a number of 
applications. When the characteristic network is used as the co-ordinate frame, the 
governing equations reduce to two ordinary differentia! equations, originally dis- 
covered by KOrrTerR (1903). Generalizations of this theory to materials whose failure 
is governed by a nonlinear relationship between the invariants p and g have been 
discussed, amongst others, by Hitt (1950) and KINGSTON and SPENCER (1970). 

The difficulties arise when one tries to relate the kinematics of the ensuing motion 
to the current stress-state. In metal plasticity the deformation is determined through 
the normal flow rule, which states that the plastic strain-rate components are pro- 
portional to the derivatives of the yield function with respect to the corresponding 
stress components. However when applied to a material yielding according to Coul- 
omb’s criterion, the normal flow rule predicts strain-rates which have a dilatational 
component which is far larger than any actually observed in the deformation of soils 
or sands. The Coulomb yield function cannot therefore serve as a plastic potential. 
HiLt (1950) noted that if the material is assumed to be incompressible, so that the 
deformation has no dilatational component, but that the principal axes of stress and 
strain-rate still coincide, then the velocity characteristics make angles of +7/4 with 
the principal directions. For such a frictional material therefore the stress and velocity 
characteristics do not coincide. This is contrary to the reasonable expectation that the 
material should fail by shearing along the lines on which the failure criterion (1.1) is 
attained. 

Two classes of theories have been developed in an attempt to overcome these 
difficulties. 


1.1. Theories of type I 


In this category of theory the Coulomb failure criterion is retained, but the normal 
flow rule is abandoned and instead the associated kinematics are developed by 
assuming that the deformation consists of a shearing motion along one family of 
stress characteristics, or more generally is the sum of two shearing motions, one along 
each of the families of stress characteristics. Theories of this general type have been 
proposed and developed by Gentev (1958), SPENCER (1964), MANDEL (1966) and DE 
JOSSELIN DE JONG (1971) and others. In these theories the material is assumed to be 
rigid/perfectly plastic and incompressible. Generalizations which incorporate dila- 
tations have been made by DE JOSSELIN DE JONG (1977), MEHRABADI and CowIn (1978, 
1981) and Harris (1985). TEUNISSEN and VERMEER (1988) have recently developed an 
elastic/perfectly plastic generalization of these models. Comprehensive reviews of 
these single and double shearing constitutive models have been written by SPENCER 
(1981) and Houtsspy and Wrortu (1982), where full lists of the relevant references 
can be found. 
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A common feature of these models is that whilst the stress and velocity charac- 
teristics are forced to coincide, the principal axes of stress and strain-rate in general 
do not. The material is nevertheless isotropic because, as emphasized by SPENCER 
(1964, 1981), the constitutive equation is actually a relation between three tensors: 
stress, strain-rate and stress-rate. For such a material it is not necessary for the 
principal axes of the first two tensors to coincide for the material to be isotropic. 

One unsatisfactory aspect of these theories is that at the present time there is no 
obvious criterion for preferring one of these various models to any of the others of 
this type. Another aspect of concern is that, as shown by SPENCER (1986), many of 
these models predict that steady shear flows are inherently unstable. 


1.2. Theories of type II 


In the second class of models, the assumption of a perfectly plastic material based 
on Coulomb’s failure criterion is abandoned and replaced by a hardening model, 
described by a family of yield curves, which depend on the current density p (or more ~ 
generally on the specific volume v). The plastic strain-rates are given by the normal 
flow rule as in metal plasticity, but the new family of yield functions are taken as the 
plastic potentials—not the Coulomb failure condition. In general as a material element 
deforms its density will increase or decrease and the local yield function will alter 
correspondingly. For plane strain deformations this family of yield curves can be 
represented by a single ‘state boundary surface’ in (g, p, p) or (q, p, v) space. 

The most successful of these models is that developed at the Cambridge University 
Engineering Department by ROSCOE, SCHOFIELD and WROTH (1958) and others and 
described by SCHOFIELD and WROTH (1968), ATKINSON and BRANSBY (1978), BOLTON 
(1979) and ATKINSON (1981). Not all the ideas of the classical Coulomb theory are 
abandoned. A key feature of this model is that the state boundary surface contains a 
so-called critical state line on which g = psinw, (or equivalently |t| = otanw,) and 
on which the direction of the normal to each of the g—p yield curves is parallel to 
the g-axis. The normal flow rule hence predicts that on the critical state line the 
associated plastic deformation occurs at constant volume. This theory is discussed 
further in the next section. 

Theories of this type, being incremental in nature, will in general require the use of 
computational packages, e.g. using finite elements procedures, to solve boundary 
value problems of practical interest. Such procedures are described by BRITTO and 
GUNN (1987). Nevertheless there is appreciable interest in developing the plane strain 
theory along the lines of classical plasticity theory for such materials in order to 
provide fundamental understanding of the nature of the deformations and hopefully 
to establish procedures for finding analytic and semi-analytic solutions in special 
cases. Thus the books by SCHOFIELD and WROTH (1968) and ATKINSON (1981) both 
contain chapters which attempt to apply Sokolovskii’s methods to boundary value 
problems, particularly in cases where the deforming material has reached the critical 
State. 

Since the normal flow rule is assumed the principal axes of stress and strain-rate 
now do coincide in contrast to theories of type I. One of the criticisms made of critical 
state models however is that as a result of this coincidence of principal axes the stress 
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and velocity characteristics will not be coincident (SPENCER, 1981). At first sight this 
would seem to be true. For example, in a region in which the material elements are 
all in a critical state, g and p are related by (1.2), but with c = 0, so that the stress 
field can be constructed by Sokolovskii’s procedures for a cohesionless soil, with the 
stress characteristics inclined at angles + (z/4—wy,/2) to the major principal direction. 
In contrast since the associated velocity field is volume-preserving it follows, as 
discussed above, that the velocity characteristics are inclined at +7/4 to principal 
directions. 

However this argument fails to recognize the incremental nature of the deformation 
of the material and that a critical state will in general only be approached asymp- 
totically. The initial purpose of the present paper is to demonstrate that the above 
argument is erroneous, that in fact when the density (specific volume) changes are 
properly allowed for the velocity and stress characteristics always coincide. In particular 
in a critical state they are both inclined at angles of +72/4 to the common principal 
directions. Essentially the same result has already been demonstrated by JACKSON 
(1983) and PITMAN and SCHAEFFER (1987), though with reference to the closely 
related problems of the dynamic flow of granular materials through hoppers. The 
methodology used here however is very different and in addition we establish the 
differential relations along the various characteristics with the aim of developing 
techniques for solving boundary value problems along the lines of classical slipline 
theory. The structure of the theory is found to be very similar to that for hardening 
models as analysed by COLLINS (1978, 1981). In addition it will be shown that there 
exists another family of weak discontinuities which are embedded in the material and 
which involve discontinuities in the derivatives of the density (specific volume) as well 
as of stress. In the critical state limit these weak discontinuities become the Sokolovskii 
characteristics. 


2. THE CONSTITUTIVE MODEL 


A soil is an aggregate of solid grains with fluid, which may be air, water or a mixture 
of both, filling the pore spaces. Soils are hence essentially two-phase materials although 
most models attempt to treat them as essentially single-phase allowing for the presence 
of the fluid by including the single pore pressure parameter. The total stress oj; is 
regarded as the sum of the effective stress o;; acting on the solid granules together 
with the pore pressure u. In a dry sand the pore pressure is zero. Boundary and initial 
value problems for wet sands and soils should properly be solved by analysing the 
coupling between the elastic/plastic deformation of the solid granules and the diffusion 
of the fluid through the porous solid matrix. Due to the difficulties of performing such 
coupled analyses attention is frequently restricted to the two limiting cases of ‘drained’ 
and ‘undrained’ loading paths. In the former the loads are applied sufficiently slowly 
for the water to drain through the matrix keeping the pore pressure constant, whilst 
under undrained conditions the loads are applied sufficiently rapidly, so that the water 
does not have time to diffuse and the resulting deformations occur at constant volume. 
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Again because of the two-phase nature of the material it is customary to work with 
the specific volume v of a material element, rather than with the density p. The specific 
volume is defined by 




























v = AV/AV, (2.1) 


where AV and AV, are, respectively, the total volume and the volume occupied by the 
solid grains in a given soil element. The value of v is hence always greater than or 
equal to 1, e = v—1 being the voids ratio. For a dry granular medium the density is 
then given by 


p = p/v (2.2) 


where p, is the density of the solid granules (assumed constant), whilst for a fully 
saturated soil, the density is 


p = (p,+(1—v)p,)/v (2.3) 


where p,, is the density of water. 
In a dry material the mass of a material element is conserved so that the density 
must satisfy the continuity or mass-conservation equation for compressible media : 


1 Dp ; 
— = 0, 4 
> Dt +divv = 0 (2.4) 
where v is the velocity vector and 
D 06 
=? = +v* grad (2.5) 


denotes the material derivative. For such a material the density and specific volume 
are related through (2.2), so that this condition can also be written 


~ — = divv. (2.6) 


On the other hand in a wet soil, the total mass of a material element is not conserved 
due to the draining away of some of the water, so that (2.4) does not apply. However 
(2.6) is still valid, since dv/v = d(AV)/AV is the volumetric strain increment in an 
element, and divv = —e, is the volumetric strain-rate (N.B. compressive strains are 
taken as positive). In developing the general equations we therefore prefer to work 
with the specific volume and (2.6) rather than use the density and the conservation 
of mass equation. 

The basic concept of Critical State Soil Mechanics is that of a state boundary 
surface, which for plane strain deformations can be represented by a convex surface 
in p, g, v space, with equation F(p, g, v) = 0 say (see Fig. la). Material elements whose 
(p,q, Vv) states lie below this surface will deform elastically in an elastic/plastic model 
or rigidly in a rigid/plastic model. These two models were termed respectively ‘Cam- 
clay’ and ‘Granta-gravel’ by SCHOFIELD and WROTH (1968). Since the analysis pre- 
sented here presumes rigid/plastic behaviour, we are using the latter class of model. 
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CRITICAL STATE LINE 








(b) O 


Fic. |. (a) State boundary surface in (p,q, v) space. (b) Family of yield loci in (p,q) space. 


The equation of the state boundary surface for ‘Granta-gravel’ is conventionally 
written 


v—T 
F(p.4.9) = yg, +inp+ ——1=0 (2.7) 


and involves three material parameters M, I and 4. This surface can also be thought 
of as generated by a family of curves in the (p,q) plane, each curve corresponding to 
a particular value of the specific volume v, the curves enlarging as v decreases (Fig. 
1b). The slope of these curves is given by 


0q/0p = q/p—M, (2.8) 


where this derivative has been evaluated at constant v. This slope is zero on the critical 
state line: 


g = Mp = sinw,p, (2.9) 


where W, = sin-'(M) is the critical state angle of friction. The variation of the 
specific volume along the critical state line is given by 


v=T—Alnp, (2.10) 


using (2.7), (2.9). 
When the state parameters of a material element reach the state boundary surface 
it will start to deform plastically and in such a way that the strain-rate components 
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satisfy the normal flow rule, with F(p,q, v) = constant as the current yield function. 
Thus on a critical state line the strain-rates have a zero dilatational component, but 
a material element compresses or expands according as its state point is to the right 
(6q/ep < 0) or to the left (6q/ép > 0) of the critical state line. Thus as the element is 
further loaded a state point on the right hand ‘wet’ side of critical moves outward to 
a yield curve with a smaller value of v, whilst the opposite occurs if the state point is 
on the ‘dry’ side and to the left of critical. The material element hence hardens in the 
first situation but softens in the second. Provided the dominant stress changes are in 
shear, in each case the state point approaches the critical state line, although on the 
‘dry’, softening side the deformation will tend to be unstable occurring in localized 
shear bands. The ‘wet’ behaviour is typical of normally consolidated clays and loose 
sands which compact when loaded. The ‘dry’ side behaviour corresponds to clays 
which are over-consolidated and to sands which are initially dense. There have been 
a number of modifications proposed to the above model. For example, ATKINSON and 
BRANSBY (1978) and ATKINSON (1981) replace the surface (2.7) on the dry side of 
critical by a Hvorslev surface which has a constant slope in the (p,q) diagram with 
equation 


F(p,q,v) = q—Hp—(M— A) exp ((T —v)/A) = 0, (2.11) 


where H is a further material parameter. The analysis presented in the following 
sections is developed for an arbitrary state boundary surface. 

Note that if pore pressures are present p and g in the above formulae relate to the 
effective stresses. 


3. UNSTEADY FLOWS 
3.1. The governing system of equations 


Since the deformations are presumed to be quasi-static the stress components @;; 
satisfy the equilibrium equations, which when referred to an arbitrary right-handed 
set of cartesian axes are 


(3.1) 





(3.2) - 


where (X,, X>) are the components of the body force vector per unit mass. Initially 
the state boundary surface will be taken as 


F (641502250 12,021,V) = 9. (3.3) 


This function must of course be symmetrical in ¢,, and o,,. Since it is written in terms 
of stress components rather than invariants the possibility of in plane anisotropy is 
allowed for. When the material is isotropic as in the Granta-gravel model, (3.3) will 
be replaced by 
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F(p,q,v) = 0. (3.3) 


In a wet soil, the pore pressure gradients, assumed known at each stage of the 
deformation, can be included in the body force terms in (3.1) and (3.2), in which case 


g;, is the effective stress tensor. 
If (v,,v>) are the velocity components, the (compressive) strain-rate tensor com- 
ponents are given by the normal flow rule: 


(3.4) 


(3.5) 


Iv Ov, ) y 
ee =e * bx, = F ; - (3.6) 


- 


where / is the proportionality factor. 
In addition the specific volume and velocity vector are related through (2.6) which 


in component form is 


a 
-< .s  Y (3.7) 


, —v —yv 
Ox, Ox, Ox, 


When analysing these equations it is necessary to distinguish between unsteady and 
steady flows. The former can be regarded as a sequence of boundary value problems 
in which (3.1)—(3.6) must be solved for the three stress components, two velocity 
components and /, the specific volume being a known function of position. Equation 
(3.7) can then be used to update the value of v for the next problem in the sequence. 
In a steady flow however 0v/0t = 0 in (3.7), which must now be solved concurrently 
with the first six equations, to determine v(x) together with the other six dependent 
variables. In this section we analyse the unsteady problem. 

The standard procedure for finding the characteristic directions and corresponding 
integrable relations of a hyperbolic system of equations is to write the system in matrix 
form and determine the eigenvalues and left-eigenvectors of a certain fundamental 
matrix (see for example HopkKINs, 1972). Whilst this procedure is universal and 
straightforward, although it can lead to some cumbersome algebra, it is essentially a 
“black box” approach, which does not bring out either much of the mathematical 
structure or the physical significance of the results obtained. Instead here we prefer 
to adopt the approach used in HILL (1950), and in a more refined form in HILL (1980), 
and argue directly from the given equations (3.1)—(3.7). 

There are essentially three different ways of defining the characteristic curves of a 
quasi-linear hyperbolic system of equations in two independent variables. 


(i) The characteristics are the curves for which the Cauchy problem has no unique 
solution. In other words, given the solution variables on a curve @, then it is not 
possible to determine the solution uniquely at a nearby point not on @. 

(ii) The characteristics are weak discontinuities, i.e. curves across which the normal 
derivatives of some of the dependent variables can be discontinuous. 
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(iii) The family of characteristics consists of curves which when used as co-ordinate 
lines are such that the original system of partial differential equations reduce to 
a system of ordinary differential equations. 


These three definitions are normally equivalent. However, in situations where the 
characteristics are repeated, which occur frequently in plasticity problems, it is not 
always possible to reduce the given system to one of ordinary differential equations. 
As will be seen in the next section, steady flow problems of density dependent granular 
materials provide such an example. The first two approaches therefore are slightly 
more general and we shall, in the main, adopt the Cauchy view point here, although 
explicit use of the other two definitions will be used when convenient. 


3.2. The characteristic directions 


The reference x,,x -axes are chosen to be locally tangential and normal to the 
proposed characteristic curve @. We suppose the dependent variables are all known, 
continuous and differentiable on @, so that all the 0/0x, derivatives are known. We 
thus have to investigate whether the given equations (3.1)—(3.6) are sufficient to 
determine uniquely all the 0/0x, derivatives, for then the solution can be found in the 
neighbourhood of @ using a Taylor series expansion. If this is not possible then @ is 
a characteristic curve. 

Looking firstly at the velocity equations (3.4)—(3.6) we note that since all the 0f/00;; 
derivatives are known on @, as is 0v,/0x,, it follows that if (3.4) can be solved uniquely 
for 4, the remaining two equations will give dv,/0x, and 6v,/0x, uniquely. Thus the 
condition that @ is a characteristic is that (3.4) cannot be solved uniquely for 7, so 
that we must have both 


Ov, 


= 3.8 
Ox, 0, (3.8) 


in which case both e,, and e,, can be discontinuous acoss @. From the latter 
equation it follows that the velocity characteristics are the zero-extension rate direc- 
tions. This is true for an arbitrary anisotropic, pressure and density dependent rigid/ 
plastic material ; it is essentially a kinematic result and would hold for any material 
whose flow rule is a relation between the stress and strain-rate tensors—the only effect 
of which would be to replace / in (3.4)—(3.6) by the distinct plastic potential function, 
g say. 

In order to consider the corresponding stress equations it is convenient to consider 
the state boundary surface function (3.3) in differential form, and following HILL 
(1980), differentiate f with respect to x,: 


f don Of bon, Ff On FF Ww _Y (3.9) 


00;,; OX, 063, Ox, 00,2 Ox, Ov Ox, 


Since the x ,-derivatives of the stress components are known on @ and the equilibrium 
equations (3.1)—(3.2) determine the x,-derivatives of o,, and o>, the only possibility 
is that 00,,/0x, may be indeterminate. In general this derivative can be determined 
from (3.9), the exception being when both 
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te ns ae 
E =0 and ke — +2— —+-= = 0. (3.10) 
007) 0072 OX) ); Iv oO: 

(Note: We are here assuming that 0v/0x, is continuous across @. Another very 
important situation arises if this derivative is allowed to be discontinuous as discussed 
in subsections 3.5 and 3.6.) 

Comparing (3.8) and (3.10) we see that the velocity and stress characteristics 
coincide and that this is a necessary consequence of the normal flow rule assumption. 
Using the equilibrium equations the second equation (3.10) can be re-written 


a, ao, Ar oo if if Ov 
Of ¢ im o = of ia ie : (3.11) 


en) - = a 
OG r> OX, O0;2 OX, 


Remembering that in an unsteady flow v(x) is known at each stage of the deformation, 
we note that the right-hand side of this equation is known, so that (3.11) is an ordinary 
differential relation for the stress components o,, and o,, differentiated in the local 
direction of @. 

If we now specialize to an isotropic material, whose state function is given in the 
form (3.3)’, the characteristic directions are given by 


Of ce l oF a (0;; —02>) OF — 


(), 3.12 
fale) il 2 Op 4q 0g { 


using (1.3). If y is the angle of inclination, measured anti-clockwise, of the major 
principal stress direction from the x,-direction, so that 


O;; =pt+qceos2y, 02> = p—qcos2y,| 


6,» = qsin 2y, j 


then it follows from (3.12) that 


011-922 


cos 2) = (3.14) 


2q 


where use has been made of the implicit function theorem. Hence the characteristic 
directions are only real as long as 


|6q/ep| = <1. (3.15) 


This condition was pointed out by Hit (1950) for a nonlinear Coulomb material, 
who showed that when it is violated the Mohrs’ circles in the (¢,t) plane lie inside 
and do not touch the yield envelope in that plane. When the characteristics are real, 
we can define a local angle of internal friction w: 


sin w = 0q/ép (3.16) 


so that y = +(n/4—wW/2). The «-/f-characteristic directions are hence inclined at 
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Fic. 2. Coincident stress and velocity characteristics, |y| = (x/4—ww/2) where sin = 0q/dp. 


angies of +(z/4—w/2) to the major principal stress directions respectively (Fig. 2). 
Whilst this is formally similar to the classical result embodied in Sokolovskii’s theory, 
it must be appreciated that in evaluating w from (3.16), the specific volume v must be 
kept constant when evaluating 0q/0p. On the critical state line g = psinyw, (2.9), but 
since v varies along this line, as given by (2.10), 6q/ép 4 siny,. In fact in a critical 
state 0g/dp = 0, so w = 0, and the coincident stress and velocity characteristics bisect 
the principal directions, in marked contrast to classical Sokolovskii theory. Since w 
is positive/negative on the “dry’’/“‘wet”’ side of critical, the characteristic directions 
are closer to the major/minor principal stress directions respectively. In the Granta- 
gravel model yw is given by, using (2.8): 


: q q , 
siny =—-—M=-—-—siny,, (3.17) 
Pp P v 


so that the characteristics are always real on the “‘wet” side of critical, but only if 
q/p < (1+sinwy,) on the “dry” side. 

Note that since y = + (2/4—w/2) the normal and tangential components of traction 
acting on elements of the a- and f-characteristics are using (3.13): 


o=p-—qsinyw and t= +qcosy. (3.18) 


3.3. Characteristic relations for velocity components 


The characteristics will in general be curved and non-orthogonal. The anticlockwise 
inclination of the major (more compressive) principal stress direction to a reference 
x-axis will be denoted by ¢ and that of the two families of characteristics by 0, and 
@, respectively, so that 


0,, 05 = $F (n/4—W/2). (3.19) 


Of course w as well as ¢# will in general vary with position. 

If (v,,v,) denote the velocity components tangential and normal to a curve, whose 
inclination to the x-axis is denoted by 6, then the condition (3.8) that the extension 
rate is zero in the tangential direction can be written 

dv, oO Ov, 


= — (y,cos@—v, sin8)|p_o = 


(3.20) 


where s measures arc length along the curve. If (v,,v,) are the orthogonal projections 
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0 
nm/2—y 


Fic. 3. Velocity projections into «-, B-curves. 


of the velocity vector in the a- and f-directionst it follows that the tangential and 
normal components on an «@-line are (see Fig. 3) 


vy, =v, and v, = vgsecy—v, tany (3.21) 


and on a f-line are 
v, =v, and v, =v, tany—v, secy, (3.22) 


so that the zero extension-rate conditions on the two families of characteristics can 
be written 


+ 


Ov, . sore 
cos as, — (vg —v, sin) _* 0, 


00, _ 
i} 





Ov 
cos W £+ (v, — vg sin) 
OS, 


When y = 0, these equations reduce to Geiringer’s equations and when yw is non-zero, 
but constant, they are formally identical with those derived by SHIELD (1953) for a 
Coulomb material with associated flow rule. Using (3.12) these equations could be 
rewritten in terms of the angular variation of ¢ and w rather than of 0, and 0,. These 
equations still apply for a material with a non-normal flow rule with a plastic potential, 
provided yw is now interpreted as the local dilatation angle—see, for example, JAMES 
et al. (1972) and FRYDMAN and BAKER (1982). 


3.4. Characteristic relations for stress components 


For general anisotropic materials the differential relation (3.11) between the stress 
components could be developed into canonical form using the general procedure 
described by HILL (1980). Here however we will restrict attention to isotropic 
materials, so that the derivatives of the state boundary function occurring in (3.11) 
can be rewritten, using (3.12)—(3.14): 


+ These are to be distinguished from the components in the non-orthogonal decomposition of v in the 
a- and f-directions. 
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If the orthogonal x,, x, axes are arbitrary and not necessarily tangential and normal 
to a characteristic, the cartesian stress components are still given by (3.13) but y is 
replaced by @. The stress gradients along the characteristics occurring in (3.11) can 
hence be written in terms of derivatives with respect to arc length as 


0o,;, OO 
Fe = py PT4C08 2M lo 


er 
- a, (4 sin 29)|4-, 


an op 
= sin 2y As +2qcos2y——, 
so that (3.11) simplifies to 


0q Ov :; 
=-2 = + p(sin 2yX, —cos 2yX,). 


The stress characteristic relations along the «- and f-curves are hence 


a 


é dq 0 
Sa 7 (3.28) 





OSp 


where the body force terms have been written in the most compact form by using 
(X,, X,) the non-orthogonal components of X when decomposed in the «-, f-direc- 
tions. In terms of the slope of the characteristics, these can be rewritten 


(3.29) 
O65 _ 
OSp “i 





When there is no density dependence and yw is constant, these relations reduce 
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formally to K6tter’s equations, and when y is variable they reduce to those given by 
KINGSTON and SPENCER (1970) for a nonlinear Coulomb material. The new feature, 
introduced by the dependence on specific volume, is the presence of known derivatives 
of the specific volume in directions perpendicular to the characteristics. This structure 
is typical of nonhomogeneous plane strain plasticity problems and may be compared 
with that of the stress equilibrium equations for an isotropic hardening material 
referred to the orthogonal sliplines which are, at an arbitrary stage of an unsteady 
deformation given by: 


0d 

OS, OSs 
A ny 7) 

c ¢ ( 
APO 


2 = 0, 
OS, 


OS, «OS, j 





using the present sign convention and notation—see COLLINS (1981). It is to be noted 
that the normal derivatives of g in (3.28) and (3.29) are not total derivatives, but only 
that part of the derivative which arises from the dependence of g on v. Some of the 
terms in these equations can be simplified somewhat for a particular state boundary 
surface such as (2.7), but their basic structure cannot be further reduced. 

These equations simplify however in the limiting critical state, for then wy = 0, the 
characteristics are orthogonal and gq and v can both be expressed in terms of p, e.g. 
by using (2.9) and (2.10). For example, if body force terms are neglected (3.29) 
become 


~ 


4) eC 


0 00 
ds. (In p)—2siny, ds, = sin, a, (In p), 


0 0 0 
as, (Inp)+2siny, ae sin w. as, (In p), 


Sp J 





where 0 = 0, = 0;—7/2. These are of course still partial differential equations and 
critical state stress fields are more readily constructed using the Sokolovskii stress 
characteristics which have associated ordinary differential relations. The relation 
between these two sets of stress characteristics is investigated in the following two 
subsections. 


3.5. Discontinuities in stress gradients 


Continuing with the notation used above in Section 3.2 we note that since equi- 
librium requires the normal (x>-) derivatives of ¢,, and o,, to be continuous across 
@, the jumps in the normal derivatives of the stress invariants p and gq are related to 
that in the derivative of o,, by 


Op 1} 00), dq . . 
45 | . 
fa | - EE 2 x |" (3.32) 


using (1.3) and (3.14), where [ ] denotes the jump in the argument. There can be no 
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such jump in the normal derivative of the third state variable v across a stress 
characteristic however, for it follows from (3.9) that any discontinuity in the normal 
derivative of v across an arbitrary curve @ must satisfy 


Of | ea), Of| ov 
ee ee (3.33) 
OG 5) OX> OV, OX) 


again since equilibrium demands that the derivatives of the other two stress com- 
ponents are continuous. Since 0f/00,, = 0 when @ is a stress characteristic (3.10), it 
follows that a discontinuity in 6v/0x, is only possible if 0f/0v = 0, which is never true 
(cf. 2.7). However it is clear from (3.33) that discontinuities in the specific volume 
gradient are possible across curves which are not stress characteristics. This possibility 
is investigated in the next subsection. 

It has been shown that the normal derivatives of p and g but not v can be dis- 
continuous across a stress characteristic. However on the critical state line, any one 
of the three state variables can be expressed uniquely as a function of each of the 
others, e.g. (2.9) and (2.10). Thus the continuity of a space derivative of any one of 
these three variables implies the continuity of the derivatives of the other two. It 
follows therefore that as the critical state is approached, the magnitude of any dis- 
continuity in stress gradient across a stress characteristic must tend to zero. In this 
sense the stress characteristics “vanish” in the critical state limit. Velocity gradients 
can still be discontinuous across these characteristics in the critical state limit however, 
so they are still “‘active” velocity characteristics. 

These results explain why in the classical direct analyses of the final failure state, 
the velocity characteristics do not show up as also being stress characteristics. Instead 
one finds that the stress gradients can be discontinuous across curves belonging to the 
completely separate family of Sokolovskii characteristics. It is shown in the next 
subsection that these characteristics are the limit of another family of weak dis- 
continuities, quite distinct from the coincident stress and velocity characteristics so 
far discussed. 


3.6. Discontinuities in specific volume gradients 


It was shown above that the normal derivative of v is always continuous across a 
stress (velocity) characteristic, but that it could be discontinuous across a distinct 
curve, providing it was accompanied by a discontinuity in the stress gradient—the 
two discontinuities being related by (3.33). However the evolutionary equation (3.7) 
also puts limitations on possible discontinuities in 0v/Ox,. Since @ is not a velocity 
characteristic, 0v,/Cx , must be continuous across this curve, so that only the first and 
third terms in (3.7) admit discontinuities, giving 


Dy Ov ov 
oa +05 = 3.34 
| | " | , 


and there is no discontinuity in divy across @. This is the condition that as time 
progresses the specific volume v continues to be continuous across a curve @ which is 
embedded in the material. A point on such a material curve will have a velocity 
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component v, in the direction locally normal to @. A general discussion of such 
discontinuity relations has been given by HILL (1961). 

The discontinuity relation (3.33) can be simplified if the material is assumed to be 
isotropic, for then 0f/da,, can be expressed as in (3.12). If 7 denotes the anticlockwise 
inclination of the major principal stress trajectory to @, (3.33) can be rewritten : 


év lat OF OF | OF 
—=—— ; 3. 
Fal ES ‘e eatin ot “ay ie 


where we have used (3.13) and (3.14), but with y replaced by 7. The relation between 
the jumps in the derivatives of o,, and p given in the first equation (3.32) is valid for 
any curve @, whilst the right hand side of (3.35) can be simplified using the implicit 
function theorem, so that finally we can write the relation between the two dis- 
continuities in normal derivatives in invariant form: 


Ov Op Ov dv | 


Thus in addition to the coincident stress and velocity characteristics there is a 
second family of weak discontinuities, which move with the material and which 
admit concurrent stress and specific volume gradient discontinuities related by (3.36). 
Although the components of the strain-rate tensor must be continuous across such a 
discontinuity, their spatial gradients need not be. This follows from the flow rule 
(3.4)—(3.6), for on differentiation with respect to x», we see that any discontinuity in 
stress-gradient will be accompanied by a corresponding jump in strain-rate gradient. 
The velocity field will hence not be continued analytically across the discontinuity. 

Now in the critical state limit we must have 


Ov Op )\dv 
— = — }— 3.37 

On On 1h ( 
on either side of the discontinuity, since now v can be expressed solely as a function 
of p, e.g. (2.10). The ratio of the discontinuities in (3.36) is hence simply dv/dp, thus 


(3.38) 


However these derivatives are also related through the identity 
dv ov a. dq ov 
dp op dp dq 


so that 


d 
4 _ siny, (3.40) 
dp 


cos 2y = 


and y = +(2/4—w,/2), which of course are the inclinations of the Sokolovskii charac- 
teristics to the major principal stress direction. We have thus established the highly 
illuminating conclusion that the Sokolovskii characteristics in the limiting critical state 
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Fic. 4. Boundary between deforming and rigid regions. 


field are the limits of weak discontinuities embedded in the material, which involve 
related jumps in the derivatives of stress and of specific volume. 


3.7. Conditions at a “‘rigid—plastic’”’ interface 


The above arguments all presume that the curve @ is embedded in deforming 
material. To complete the picture we must investigate the situation where the material 
is rigid on one side of this curve, so that @ is part of the boundary of the deforming 
zone. As before we choose axes Ox,, Ox, locally tangential and normal to @, with 
Ox, directed into the deforming zone, so that f= 0 on this side of @. At a point on 
the rigid side, however, the state (p,qg,v) point may lie inside the state boundary 
surface and so we only have f < 0 (Fig. 4). 

Since in general? the strain-rate is discontinuous across such a bounding curve, @ 
must be a velocity characteristic and so 6f/00,, = 0 (3.8). In order to discuss the local 
stress field we again use the equilibrium equations (3.1) and (3.2), which of course 
hold on both sides of @, and the differentiated form of the state boundary surface 
condition (3.9), except that now on the rigid side of @, the equality in (3.9) must be 
replaced by the inequality >0. Writing this slightly modified form of (3.9) for both 
sides of @ and subtracting gives the conditions for the existence of discontinuities. 
Since 6f/0a,, = 0 and the equilibrium equations demand that the second and third 
terms are continuous across @, we conclude that 


Ov 
<0 (3.41) 


Ox deforming Ox rigid 


provided 0f/dv > 0, as in the Granta-gravel model. 

Thus both the normal derivatives of o,; (or equivalently p) and vy can be dis- 
continuous at a rigid/deforming interface, but the jumps are not related as in (3.33). 
Instead the jump in the specific volume gradient has to satisfy the inequality (3.41). 


+ There is in fact another possibility. If the strain-rate is continuous, but its gradient is discontinuous 
across @, then the bounding curve could be a member of the second family of embedded weak discon- 
tinuities. In the critical state limit, such a bounding curve would be a Sokolovskii stress characteristic. 
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4. Strapy FLows 


Due to the problem of having to solve simultaneously for the motion of the pore 
water it is really only sensible to discuss the steady flows of dry granular materials in 
the present context. For such materials the specific volume v can be replaced by the 
density p in the evolutionary equation, cf. (2.4) and (2.6). However to be consistent 
with the rest of the paper we will continue to use the specific volume as the state 
parameter. 


4.1. The characteristic directions 


The analysis of the stress equations (3.1)—(3.3) and velocity equations (3.4)—(3.6) is 
exactly the same for steady and unsteady flows. The stress and velocity characteristics 
coincide and are inclined at angles of +(/4—w/2) to the major principal stress 
trajectory. However we now must simultaneously solve the evolutionary (continuity) 
equation (2.6) which for steady flows reduces to: 


(4.1) 


We proceed to investigate the characteristic structure of this equation using the 
same procedure as for the others, by assuming that the x,-axis is locally tangential to 
a characteristic curve @. We note that provided © is not a velocity characteristic, the 
normal flow rule equations (3.4) and (3.5) enable the derivative @v,/0x, to be written 
in terms @v,/0x, which is in the direction of @, viz. 


=. [a= = ). (4.2) 
CX, CO >? CG) 


Thus if v, = 0 on @, the remaining three terms in (4.1) form an ordinary differential 
relation along this curve. In other words, the streamlines form a fifth family of 
characteristics. This result is otherwise obvious, since for a steady flow the streamlines 
coincide with the particle paths and the evolutionary or continuity equations (2.6) or 
(2.4) are just expressions for the material derivatives of specific volume or density, 
respectively. Note however that although (4.1) would allow a discontinuity in the 
normal derivative 0v/Ox, across a streamline, the differential form of the state bound- 
ary yield condition (3.9) shows that this is not possible unless 0f/év = 0. Since this 
condition is never true in the critical state models, it follows that any discontinuity in 
v across a streamline must be in a derivative of higher order than one. Further analysis 
Shows that none of these conclusions is altered in the exceptional circumstance of a 
streamline coinciding with a velocity characteristic. 


4.2. The characteristic relations 


Putting v, = 0 in (4.1) and substituting for 0v,/0x, from (4.2) the ordinary differ- 
ential relation along the streamline can be written: 














Characteristics theory for soils 








Fic. 5. Orientation of streamlines. 


Ov = —y¥ = = 0 (4.3) 
Os Os 
or 
O O 
’ Os a Os (ing) = 0, (4.3)’ 


where v( = v,) is the speed, 0/ds = 0/0x, denotes differentiation along the streamline 
(note that since v, = v, = 0 there is no curvature term as in (3.20)), and 6 is the ratio 


0 = Of/00,,/(0f/00,, + 0f/0022) = e);/(€1; +e22). (4.4) 


These state boundary function derivatives can be simplified when the material is 
isotropic as in (3.12), in which case using (3.13)—(3.15), 6 can be written 


.  sinyy—cos2y _ sing, SIN Xp 


2siny ——-cos(y,+X,)’ 9) 


where x is the anticlockwise inclination of the major principal stress trajectory from 
the streamline direction, and y,, 7%, are the angles between the streamline directions 
and a- and f-sliplines respectively (Fig. 5). Note that 6 = 0 when the streamline 
coincides with a slipline, so that the speed is then constant along the streamline, 
which of course is necessary since a slipline is a zero-extension rate direction. If the 
streamlines are inclined at 2/4 to the major principal stress trajectory, 6 = 1/2 and 
v*/y is invariant along a streamline. In general 6 will take a constant value and v/v? 
will be invariant on a streamline in flows in which the product tan x, tan xy, remains 
constant. Physically this occurs when the ratio of the extensional strain-rate in the 
streamline direction to the dilatation rate is a constant. This is not the case in the 
important special class of diagonal streaming flows (HILL, 1966) in which y = 0 and 
6 = —cot(/2). 

The characteristic relations along the coincident stress and velocity characteristics 
are again given by (3.28) or (3.29) and (3.23). However the former involve derivatives 
of the specific volume v in the direction normal to the characteristics. In unsteady 
problems these derivatives are known at each stage of the deformation, but in a steady 
flow we have to solve for v as one of the dependent variables. Equations (3.28) and 
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(3.29) are hence, in the case of steady flows, not ordinary differential equations as 
one normally expects on characteristics. This situation is typical of steady flows of 
hardening materials (see COLLINS, 1981), and is a consequence of the system of 
equations having repeated characteristics. This can perhaps be best understood in 
terms of the standard matrix formulation of the characteristics theory of quasi-linear 
systems of equations, and is discussed here in the Appendix. It should be remembered 
that steady flows of hardening materials are not well-posed problems; in order to 
guarantee uniqueness the deformation must be followed incrementally (HILL, 1956). 


5. DISCUSSION 


The analysis in the previous sections has established that two types of “weak 
discontinuity” can arise in the quasi-static, plane strain deformation of an isotropic 
rigid/plastic material, whose yield condition depends on the specific volume and which 
deforms according to the normal flow rule. 

(a) Coincident velocity and stress characteristics, which are inclined at angles of 
+(z/4—w/2) to the major principal stress directions, where sinw = 0q/dp. (In an 
earlier attempt at such an analysis by JENIKE and SHIELD, 1959, the principal axes of 
stress and strain-rate were assumed to coincide, but since the deformation was assumed 
to be always effectively isochoric, the normal flow rule did not apply so that these 
authors found that the stress and velocity characteristics did not coincide.) The 
gradients of the velocity and stress invariants, but not the specific volume can be 
discontinuous across these characteristics. The characteristics are in the zero-exten- 
sion-rate direction and hence admit the ordinary differential relations (3.23). The 
corresponding differential relations for the stress invariants (3.28) or (3.29) involve 
the derivatives of the specific volume in the direction normal to the characteristics. 
These latter derivatives are known at a generic stage of an unsteady deformation, but 
have to be found as part of the solution in a steady flow. In the critical state limit, 
Ww = 0 and these characteristics bisect the principal directions and coincide with the 
maximum shear stress (strain-rate) trajectories. In this limit the velocity characteristics 
remain “active”, but the magnitude of all discontinuities in the stress gradients tends 
to zero and these curves are no longer “active” stress characteristics. 

(b) Discontinuities in the normal derivatives of the specific volume can occur across 
curves which are embedded in the material. Such discontinuities must be accompanied 
by discontinuities in the pressure gradient as prescribed by (3.36). In the critical state 
limit these discontinuities become the classical stress characteristics of Sokolovskii 
theory inclined at angles of +(2/4—w,/2) to the major principal stress directions. 

If the soil or granular material is initially on the “dry” side of critical it will dilate 
and soften. It will hence become unstable and deform in narrow dilating localized 
shear bands. A body made of such material hence “‘fails’’ before the critical state is 
achieved. However if the soil is on the “wet” side of critical it will compress and 
harden, the overall distortion is small and velocity discontinuities cannot form until 
the critical state is reached. For material on the “wet”’ side, therefore, the limiting 
critical state represents the ultimate failure state. 
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FiG. 6. Sokolovskii stress (dotted lines) and velocity (full lines) characteristics for passive failure of soil at 
a vertical retaining wall. 


Due to the difficulties in computing the continued solution to a boundary value 
problem for a hardening material, there is great interest in bypassing the preliminary 
deformation and going directly to the final critical state failure solution. Thus the 
books by SCHOFIELD and WRoOTH (1968) and ATKINSON (1981) have chapters dis- 
cussing techniques for directly computing the failure loads by assuming all the material 
is already at the critical state. Possible stress fields are most easily computed using 
Sokolovskii’s methods starting from a surface on which both components of traction 
are known (e.g. a free surface). The corresponding orthogonal velocity characteristics 
which make an angle of y./2 with the Sokolovskii characteristics can then be drawn. 
A representative example which is the solution for the passive failure of a retaining 
wall is shown in Fig. 6 (cf. BRANSBY, 1973). The constant state stress region ACD is 
continued through the centred fan ACB to the wall. When the wall rotates about B, 
a discontinuity in the velocity gradient is generated along the velocity characteristic 
BEF, which represents the boundary of the deforming region across which the stress 
and specific volume gradients may be discontinuous, but will satisfy (3.41). The 
stresses in the rigid region outside this curve are not unique, they may not be at yield 
and so need not be compatible with the Sokolovskii field. The line AE is a discontinuity 
in the gradients of stress, specific volume and strain-rate and is an example of the 
limit of the type of embedded discontinuity discussed in Section 3.6. The presence of 
the region between the outer velocity and Sokolovskii stress characteristics which has 
caused much concern and discussion in the past (e.g. FRYDMAN and BAKER, 1982) is 
seen to be of no particular physical significance. The velocity characteristic BEF is 
the actual boundary of the deforming region. It is convenient to construct the complete 
Sokolovskii network initially, but once the position of the boundary of the deforming 
zone has been determined, the stress field outside this boundary can be disregarded. 

Once the stress field has been determined using the Sokolovskii characteristics, the 
specific volume v can be calculated at each point in the deforming zone using (2.10). 
The asymptotic distribution of v in this zone at failure is hence, within the limitations 
of small strain theory, independent of its initial distribution or loading history. 

Critical state stress solutions are also of interest from the viewpoint of the appli- 
cation of the extremum principles of limit analysis. Since all the state points on the 
“wet” side of critical lie below the inclined plane, with equation 


q—psiny, =0 (5.1) 


in (p,q,v) space, it follows that their stress states are statically admissible for a 
perfectly plastic material with (5.1) as its yield condition. It follows therefore that in 
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any boundary value problem the failure load for such a perfectly plastic material must 
be greater than, or equal to, any of the loads which correspond to the material being 
anywhere on the “‘wet”’ side of critical or at the critical state. Hence such an idealized 
solution provides an upper bound to all the successive yield loads, including the 
ultimate failure load, providing this is approached from the ‘“‘wet”’ side of critical. 
However this limit theorem can only be applied if the velocity field associated with 
this perfectly plastic solution is given by the normal flow rule with the yield condition 
(5.1) as the plastic potential. In this context therefore the associated velocity field will 
dilate, with y, as the dilatation angle, and the velocity characteristics now coincide 
with Sokolovskii’s stress characteristics. Thus for bounding purposes the theory 
developed by SHIELD (1953) is appropriate but when the critical state is viewed as an 
asymptotic state the isochoric velocity field theory described by HILL (1950) is the 
correct model to choose. It should also be noted that if the complete solution for the 
idealized perfectly plastic problem cannot be found, approximate kinematic upper 
bound solutions will still give over-estimates to the actual failure load. However 
solutions of lower bound type, using statically admissible stress fields,will not neces- 
sarily provide under-estimates of the actual failure load. For this reason it can be 
argued that upper bound approaches are more appropriate—see, for example, ATKIN- 
SON (1981) and COLLINS ef a/. (1988). An important exception, however, is the 
undrained loading of normally consolidated clays, where the specific volume and 
hence the effective yield stresses do not change, so that the material can be effectively 
treated as a metal. 

Finally, although the fact that the stress and velocity characteristics do after all 
coincide for density dependent hardening models, this fact does in no way disprove 
the validity of any of the competing perfectly plastic models—labelled type I in the 
Introduction—however it does perhaps remove one of the original motivations for 
developing these theories. The major remaining difference between the two classes of 
theories is that the perfectly plastic models require non-associated flow rules, so that 
the principles axes of stress and strain-rate do not coincide, whilst the critical state 
family of models retain the normal flow rule assumption so that the principal axes do 
coincide. There appears to be no clear experimental evidence which would settle this 
issue conclusively. It must be borne in mind that both families of models assume 
the material is isotropic whereas real granular materials soon develop pronounced 
anisotropic behaviour when deformed. In addition it does not seem possible to predict 
all types of observed behaviour, such as shear band localization, without complicating 
the models further, e.g. by the introduction of yield vertices. 


6. CONCLUSION 


Many of the difficulties and paradoxes which have arisen in trying to develop the 
kinematics of plane strain deformation for granular materials are due to the desire to 
make the velocity characteristics coincide with the stress characteristics, which in the 
past has always meant the Sokolovskii stress characteristics. The analysis presented 
here shows that in the critical state soil model this paradox is resolved in a natural 
and elegant way. There are in fact two families of weak discontinuities, one having 





Characteristics theory for soils 23 


the orthogonal velocity characteristics as their critical state limit, the others yielding 
the Sokolovskii stress characteristics in this limit. These two families have a different 
physical significance in the deformation prior to the critical state and there is no 
necessity for them to coincide in the final failure state. 
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APPENDIX 


A quasi-linear system of equations in two independent variables can be conveniently written : 


Au,+ Bu, = f, (Al) 


where u is the n-dimensional column vector of solution variables, A and B are (n x n) matrices 
and f is a given column vector. The elements of A, B and f may depend on u, x and y but not 
on the derivatives of u. We can form a linear combination of these equations by premultiplying 
by a row vector ¢’. The resulting single equation will be an ordinary differential relation in the 
direction which makes an angle @ with the x-axis if there exists a second row vector m’ such 
that 
¢'(Au,+Bu,) = 7'f 
= mu, = m’(u, cos0+u, sin @) 
so that 
f’A=cosOm’ and /’B=sin0m’. 


Eliminating m’ gives the eigenvalue equation 
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4'(A—/B) =0’, (AS) 


where the eigenvalues A = cot@ give the characteristic directions. The corresponding left- 
eigenvectors ¢’ give the multipliers by which the original equation has to be multiplied to give 
the ordinary differential relation along the characteristic. However if a particular eigenvalue is 
repeated twice, then the corresponding eigenvector may be unique or it may belong to a two- 
dimensional subspace (WILKINSON, 1965). In the latter situation it is not possible to find an 
ordinary differential relation along the corresponding characteristic. The repeated eigenvalues, 
giving the coincident stress and velocity characteristics, have eigenvectors which are unique in 
the case of a perfectly plastic material, but which are non-unique for the steady flows of a 
hardening material. 
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ABSTRACT 


A WELL-KNOWN class of constitutive relations for cohesionless granular materials are the double-sliding 
models. The double-sliding model proposed by SPENCER (J. Mech. Phys. Solids, 12, 337, 1964) for 
incompressible materials as well as those of SPENCER and KINGSTON (Rheol. Acta, 12, 194, 1974) and 
MEHRABAD!I and CowIn (J. Mech. Phys. Solids, 26, 269, 1978) for compressible materials are generalized 
to describe the unsteady flow of cohesionless granular materials. These models are analysed and it is shown 
that they lead to equations that are not well-posed. Implications of this finding are discussed. 


1. INTRODUCTION 


THE STRESSES in cohesionless granular materials are reasonably well described by the 
Coulomb yield criterion. Various theories have been proposed to describe the kine- 
matic behaviour of granular materials. A well-known model for large, continued 
deformations in two dimensions was developed by SPENCER (1964). 

It postulates that the material is isotropic and incompressible and satisfies the 
Coulomb yield criterion. The equations that govern the kinematic behaviour are 
derived from physically plausible assumptions. The deformation is assumed to occur 
by a superposition of two simple shearing motions along the sliplines (where the 
critical ratio of shear to normal stress is attained) and a rotation of the slipline 
field. It was from these principles that SPENCER (1964) derived his so-called velocity 
equation. 

In Spencer’s original model it is assumed that the material is incompressible. Various 
generalizations have been proposed to include the effect of the compressibility of the 
material: SPENCER and KINGSTON (1974) and MEHRABADI and CowlIn (1978). In its 
simplest form the proposal of SPENCER and KINGSTON (1974) states that the density 
is determined by the mean pressure, whereas MEHRABADI and CowIN (1978) included 
the effect of the dilatancy, i.e. the volume change induced by shear deformation, into 
the double-sliding model. 

In contrast to plasticity relations based on a smooth flow potential the double- 
sliding models predict that the stress and deformation-rate tensors are non-coaxial. 
Experimental evidence (DRESCHER and DE JOSSELIN DE JONG, 1972 ; ODA and KONISHI, 
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1974: DRESCHER, 1976; ALLERSMA, 1985) confirms the existence of this non-coaxiality 


(at least under certain conditions). 
Spencer’s incompressible model has been generalized by MORRISON and RICHMOND 


(1976) to describe the unsteady flow of granular materials. The same will be done 
here for both models mentioned that include the compressibility of the material. 

The complete system of equations that will be studied here consists of the momen- 
tum equations, the continuity equation, the Coulomb yield criterion, the velocity 


equation and a density equation. 

It will be shown that these equations do not satisfy the Hadamard stability criterion 
(JOSEPH et al., 1985), which requires that plane-wave solutions remain bounded as 
t— o. This means that the double-sliding models are not well-posed. In the case of 
Spencer’s original incompressible model this was previously stated by SCHAEFFER 


(1987). 
Tensile stresses will be reckoned as positive, this being the usual convention in 


continuum mechanics. The stress tensor is denoted by T and the velocity vector by 
u = (u,v)’. The rate of deformation tensor D and the vorticity tensor W are defined 


by 
D = 3[(Vu)+(Vu)’], (1.1) 
W = 3[(Vu) —(Vu)’]. (1.2) 
In two dimensions the deviator D’ of D is given by 
D’ = D— \(tr DDI. (1.3) 


The deviator T’ of the stress tensor is similarly defined. The mean pressure p is given 
by 
p=-—itrT (1.4) 


and the invariant t of the stress deviator by 


t* =4tr (T’°T’). (1.5) 


2. THE BAsic EQUATIONS AND THE COULOMB YIELD CRITERION 


The momentum equations in the vertical and horizontal direction are respectively 


(2.1) 


) a ale 
| 


The continuity equation reads 





eo 
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0p Op Ou Ov 
+u +9 + = (). (2.3) 


Ot Ox Ox doy 


It is assumed that the material is cohesionless and satisfies the Coulomb yield 
criterion, which can be expressed in the Sokolovski variables y and p as 


E xs aaa — pl —sin (p) COS (2y)I, (2.4) 
T,, = —pll +sin () cos (2) (2.5) 
T,, = p sin () sin (2y). (2.6) 


Here y is the angle of the direction of the algebraically greater principal stress with 
the (vertical) x-axis and ¢ is the internal friction angle. 


3. THE KINEMATIC BEHAVIOUR 


The generalized double-sliding models consist of two additional equations that 
describe the kinematic behaviour: a velocity equation and a density equation. 
The generalized velocity equation as proposed by MEHRABADI and CowIn (1980) 
is 
T:-T—T-T = 2h,[T-D—D-T}. (3.1) 
In (3.1) T denotes the objective Jaumann derivative 


T =T-—W-T+T-W. (3.2) 


In (3.2) the superimposed dot denotes the material time derivative. The parameter /, 
in (3.1) is the so-called non-coaxiality parameter. It is a measure of the non-coaxiality 
between the stress and deformation-rate tensors. In scalar form (3.1) reads 


W - W., - =! [2 COS (2y )D., —sin (2p MD, eae Dy, )]. (3.3) 


Two density equations have been proposed: a simple relation between density and 
pressure was given by SPENCER and KINGSTON (1974) and a more sophisticated 
dilatancy relation by MEHRABADI and CowIn (1978). 

SPENCER and KINGSTON (1974) proposed the following relation : 


dp 
dp =F (oPoku). (3.4) 


Here k,, denotes a hardening parameter that describes the stress and deformation 
history. In its simplest form, if the effect of the stress and deformation history is 
neglected, (3.4) reads 


p=f(p). (3.5) 


It will be assumed here that the density is a strictly increasing function of the mean 
pressure p. 
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The Mehrabadi—Cowin dilatancy relation was derived by MEHRABADI and COWIN 
(1978) and Harris (1985). It gives a relation between the volume increase induced by 
shear deformation and the distortional work. This dilatancy equation is 


ttr (D) =f tr (T’-D’). (3.6) 


The dilatancy parameter f is defined by f = sin (v)/cos (@—v), where v (v < #) is the 
dilatancy angle. The original dilatant double-sliding model of MEHRABADI and COWIN 
(1978) is retrieved if /, in (3.1) is given by 


h, - 1—f sin (¢) 
t  f-—sin(¢) 


The incompressible model of SPENCER (1964) is a special case of the Mehrabadi 
Cowin model: it is obtained if vy = 0 and hence f = 0. 


(3.7) 


4. ANALYSIS OF THE GENERALIZED SPENCER—KINGSTON MODEL 
The generalization of the Spencer—Kingston double-sliding model to the unsteady 


flow of cohesionless granular materials consists of (2.1)—(2.6), (3.3) and (3.5). Sym- 
bolically these equations can be written, after the substitution of (3.5) into (2.3), as 


Ow Ow 
= A(w) — +B(w)— + 
Ox Oy 














where w = (u,v, ,p). The matrices A and B are given in Appendix A. 
A plane-wave solution will be substituted into the system of equations (4.1). This 
plane-wave solution has the form 


W = Wyo exp [i(@f+k>x)], (4.2) 


where k = (k cos (9), k sin (9))’ and x = (x, y)’. 
Disregarding the source term it follows from (4.1) that non-trivial plane-wave 
solutions exist if and only if 


det cos (9)A+sin (9)B— - ‘ (4.3) 


where I is the unit matrix. 

Equation (4.3) can be solved for w. If all solutions are real (for all 9) then the 
system of equations is hyperbolic. However, if for some 9, one of the solutions has a 
negative imaginary part, then the problem is not well-posed. In this case the solution 
can grow exponentially with time, as can be seen from (4.2). This means that the 
equations do not satisfy the Hadamard stability criterion. 
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The solutions for w are given in Appendix B, where it is also shown that no value 

of h, exists for which the Hadamard stability condition is satisfied. This means that 
the Spencer—-Kingston double-sliding models are not well-posed. 


5. ANALYSIS OF THE GENERALIZED MEHRABADI-COWIN MODEL 


- 


The generalization of the Mehrabadi-Cowin double-sliding model to the unsteady 
flow of cohesionless granular materials consists of (2.1)—(2.6), (3.3) and (3.6). Sym- 
bolically these equations can be written as 


Ow Ow 
= C(w) — +D(w)— + 
Ox oy 














— — 


where w = (u,v, /, p, p). The matrices C and D are given in Appendix C. 
Again a plane-wave solution will be substituted into the system of equations. This 
plane-wave solution has the form 


W = Wo exp [i(w@/+k->*x)], (5.2) 
where k = (k cos (3), k sin (9))’ and x = (x, y)’. 
Disregarding the source term it follows that non-trivial plane-wave solutions exist 


if and only if 


~ 


det cos (9)C+sin (9)D— 


W 

k 

where the modified unit matrix I is defined by 

=|! wins and i<5 (5.4) 
0 otherwise. 


The Hadamard stability condition again requires that the imaginary part of the 
solutions for @ is never negative. 

The solutions of (5.3) for w are given in Appendix D, where it is also shown that 
no value of /, exists for which the Hadamard stability condition is satisfied. This 
means that the dilatant double-sliding models are not well-posed. Therefore the 
Statement of MORRISON and RICHMOND (1976), that the generalization of Spencer’s 
original model to the unsteady flow of granular materials leads to a hyperbolic system 
of equations, is incorrect. This was also remarked by SCHAEFFER (1987). 
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6. DISCUSSION 


It has been shown that the incompressible and dilatant double-sliding models do 
not satisfy the Hadamard stability condition. This means that, according to these 
models, the material always shows unstable behaviour. For a simple shearing flow 
this was also shown by SPENCER (1986). 

Stable flows of granular materials are physically possible, although they may be 
difficult to produce experimentally as there is a tendency for the formation of shear 
bands. Hence it is concluded that the deformation behaviour of granular materials is 
not described satisfactorily by the double-sliding models. 

Apparently the double-sliding models are too restrictive. Possibly a better descrip- 
tion of the material behaviour is given by the double-sliding free-rotating models of 
DE JOSSELIN DE JONG (1971, 1977). In these models (3.3) is replaced by 


Q + W., ” =! [2 COS (2W)D,, — sin (2W)(D,. - D,,), (6. I ) 


where /, is given by (3.7) and Q is an undetermined rotation rate. 

Some experimental support for these models is given by DRESCHER (1976). A 
disadvantage of the free-rotating double-sliding models is the indeterminacy of the 
rotation rate Q, which makes it impossible to perform definite calculations. 

Micro-mechanical modelling (CHRISTOFFERSEN ef a/., 1981) could lead to a relation 
between Q and other quantities describing the stresses and deformation of the granular 
material. An alternative to this approach is the numerical modelling of assemblies of 
particles (CUNDALL and STRACK, 1979). In this way it may be possible to obtain a 
model satisfying the Hadamard stability condition, but lacking the indeterminacy of 
the free-rotating double-sliding model. 
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APPENDIX A 


The matrices A and B from the system of equations (4.1) are given by 

















é: baal wW)7 
—u 0 thins <= 
p p 
l 
0 = 2? sin () cos (2) sin (#) sin (2) 
i a p p 
—h, . " 1, 5 0 
> sin (2) 2+ 5 cos ( yy) —u 
ee 0 0 ~u 
= p al 
and 
q ¥ - 7 
—v 0 " sin (@) cos (2) r sin (~) sin (2y) 
—1l-si s (2 
0 —v 2? sin (d) sin (QW) — vod = (2) 
B= p p 
| h, h, ° 
—5,+ = cos (2yW) =sin (2~) —v 0 
2t 2t 
—p 
q 0 P 0 —1 : 
where p’ is dp/dp. 
APPENDIX B 
Define z = —(mw/k)—cos (9)u—sin (9)v. 
It follows after some lengthy algebra that the solutions of (4.3) are determined from 
z*—Az*+B=0, (Bi) 


where 
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h l 
A -Psin (@)| ' + cos (2 — 29) | + [1 —sin (@) cos (2 —29)} 
p t p 


h : 
B= P sin ()| 1+ ‘cos (2 28) [feos 2p -29)—sin ()}. 
pp T 


So z 1s determined from 
z* = ${A4+(A*—4B)"’). (B4) 


The solutions for z have a non-negative imaginary part if and only if: A > 0, B > 0 and 
A*—4B > 0. It is impossible to find 4, such that all solutions for z are always real or have a 
positive imaginary part. This means that the generalized Spencer—Kingston double-sliding 
model always exhibits Hadamard instabilities. 


APPENDIX C 


The matrices C and D from the system of equations (5.1) are given by 
i 


>. —1+sin (@) cos (2¥)7 
—u 0 _2! sin(d)sin(2yv) 0 p 
) p 


f 
. | 
2~ sin (@) cos (2) , Sin (p) sin (2) 
p f 


ae ee ‘ 
>, Sin(2y) 34+ cos (2p) 
T T 


—p 0 
1 —f cos (2) — B sin (2y) 








2? sin () cos (2) - sin (#) sin (2W) 
pP f 


—1—sin(@)cos (2y) 


2? sin (#) sin (2W) 
p pP 


| h, — hy : 5 
—5+ 7, 008 (2H) 5, sin (2) 


—_ ye 


0 om 
L —f sin (2p) 1+ fcos (2) 








APPENDIX D 
Define z = —(w/k)—cos (9)u—sin (9)v. 
It follows after some lengthy algebra that the solutions of (5.3) are determined from 
Cz’ = Dz, (D1) 
where 
C = [(1+£ sin (#))—(sin (¢) + B) cos (2y —29)] (D2) 


and 
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h, l 
D = 2 sin (#)[{sin (@)—cos ey—20)]| - y+ +(e 1) cos ey-29)| (D3) 


Since v < ¢ it follows that C > 0. The solutions for z have a non-negative imaginary part if 
and only if D > 0, which is only possible if 
h h ; 
-414— B+ (\p- s) cos (2y -28)| = K[sin (@)—cos (2y —29)], (D4) 


2t 


where K should be positive. It follows from (D4) that 
h, 1—fsin(@) 
t  f—sin(¢) 
Hence it is impossible to find 4, such that all solutions z are always real or have a positive 
imaginary part. This means that the generalized dilatant double-sliding model always shows 
Hadamard instabilities. 


h 
and - —p> 0. (D5) 
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ABSTRACT 


FOR CREEP crack growths by the nucleation and growth of grain boundary cavities, steady-state crack 
growth rates (da) are derived when the cavity nucleation is assumed to be continual and strain controlled. 
For the crack growth under HRR fields, 4 is linearly proportional to C* as reported in other models. 
However, when the crack growth occurs under elastic fields, d ~ Kj?**”° (n > 8/3) for diffusional cavity 
growth and d@ ~ Kj (n > 2) for cavity growth by power law creep. Predictions of the model are applied to 
a Ni-Cr steel and a Nimonic 80A, and compared with other models through analyses on the dependence 
of a on macroscopic load parameters and the activation energy of 4. 


NOMENCLATURE 


crack growth rate 

path independent integral for creeping solids 
stress intensity factor (mode I) 

stress tensor 

equivalent stress 

polar coordinates centered at the crack tip 

creep strain rate tensor 

coefficients in Norton’s creep equation (é = Bo”) 
angular function in elastic singular solution, (1) 
angular function in HRR singular solution, (3) 
dimensionless parameter in HRR solution, (3) 
number of cavities per unit grain boundary area 
time derivative 

stress exponent in the cavity nucleation rate, (4) 
temperature dependent coefficient in the cavity nucleation rate, (4) 
volumetric growth rate of a cavity 

grain boundary area occupied by a cavity 

grain boundary area occupied by a cavity at nucleation 
atomic volume 

grain boundary diffusivity 

grain boundary thickness 

Boltzmann’s constant 

temperature 

cavity tip dihedral angle 
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functions of y describing the cavity shape 

cavity radius 

dimensionless factor representing triaxial stress states in (8) 
initial crack length 

crack length increment 

cavitated area fraction at grain boundary 

critical damage for the crack advance 

damage accumulated during the incubation period 

damage accumulated during the crack growth 
characteristic distance 

incubation time for the crack extension 

critical tensile stress for the cavity nucleation 

distance from the crack tip where o,, = a, 

distance from moving crack tip 

dimensionless factor in (24) 

dimensionless factor in (36) and (37) 

parameter showing the density of cavitating facets 
additional factor to /, caused by cavitating facets in HRR singular solution in (52) 
melting temperature, K 

Young’s modulus 

activation energy of crack growth rate [a = d)exp(—Q,/RT)]| 
activation energy of grain boundary diffusion 

activation energy of creep [é = é) exp (—Q,/RT)] 


INTRODUCTION 


UNDER CREEPING conditions, crack tip deformation fields are time dependent and a 
crack usually propagates along the grain boundary by the nucleation, growth and 
inter-linkage of grain boundary cavities. Thus, the creep crack growth problem has 
been studied extensively on the basis of the singularity analysis and the micromech- 
anism of grain boundary cavitation (HUTCHINSON, 1968; RICE and ROSENGREN, 
1968; Riepet and Rice, 1980; Hur and RiepeL, 1981). In an elastic/power law 
creeping material, RmepEL and Rice (1980) showed that crack tip fields are char- 
acterized by the elastic stress intensity factor K, right after load application but 
ultimately relax to steady-state HRR fields characterized by the C* integral as creep 
strain accumulates near the crack tip. However, as long as the region of creep strain 
dominance, the so-called creep zone, is much smaller than the crack length, K, is the 
appropriate load parameter to be correlated with crack tip fields and the creep crack 
growth rate d. This is often stated as the small scale creep (SSC) condition. On 
the other extreme, under extensive creep conditions where the creep zone becomes 
comparable to the crack length or specimen dimensions, C* is the characterizing load 
parameter. Experiments show that ad ~ C* in many ferritic steels while a ~ (K,)’ 
in Ni-base superalloys with p ranging between 2 and 9 (FLOREEN, 1975, 1983: 
SADANANDA and SHAHINIAN, 1977, 1978, 1983; NAzMy and WuTuricu, 1983; 
RIEDEL, 1983; RIEDEL and WAGNER, 1984). Superalloys generally manifest a very 
long transition time between the SSC and extensive creep. 

With recent developments in our understanding of the creep cavitation process, 
micromechanistic models of creep crack growth are proposed which incorporate 
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cavity nucleation and growth under crack tip fields (WILKINSON and VITEK, 1982; 
RIEDEL, 1983; EARTHMAN, GIBELING and Nix, 1985). Often, cavities are assumed to 
nucleate at the attainment of a critical local stress and grow without further cavity 
nucleation during the subsequent crack advance. However, it is well documented in 
literature that cavities nucleate continually in the near tip region of cracked specimens 
as well as in smooth rupture specimens, because fractographs of the crack tip region 
reveal cavities of varying size with different nucleation times (YU, CHO and PARK, 
1989; Yu and HONG, 1989). Furthermore, in a study on creep rupture, the writers 
have shown that the continual cavity nucleation profoundly influences the rupture 
time, which suggests a definitive work to be done on the role played by the continual 
cavity nucleation in the creep crack growth (YU and CHUNG, 1988, 1989). 

Here, we present a steady-state creep crack growth model when the crack tip fields 
are characterized by either K, or C*. The cavity nucleation is assumed to be initiated 
at the attainment of a critical local stress and continues in a strain controlled way 
until the local accumulated damage reaches a critical value. In determining the cavity 
nucleation rates, experimentally observed linear relationships between the grain 
boundary cavity number and the creep strain are utilized (DYSON and MCLEAN, 1972; 
CANE and GREENWOOD, 1975; FLECK, TAPLIN and BEEveRS, 1975; CANE, 1979; 
NEEDHAM and GLADMAN, 1980). After nucleation, cavities are assumed to grow by 
diffusion or creep under the stress fields characterized by K; or C*, and any interactions 
among cavities during growth are ignored. With respect to the crack tip which moves 
with a constant velocity a, crack tip stress fields are assumed to be time independent 
and any additional effects caused by the evolution of cavitational strain on the 
asymptotic fields are ignored. 


2. CRACK TIP FIELDS AND CAVITATION BEHAVIOR 
2.1. Crack tip stress fields 


In Ni-base superalloys, size of the creep zone is usually much smaller than the 
specimen size and a crack grows predominantly in elastic surroundings where crack 
tip stresses are given by 


g;(r,0) = * f(0). (1) 


2ur 


Here, o;; is the stress tensor, (r,@) are polar coordinates centered at the crack tip so 
that r is the distance from the crack tip and 0 = 0 is the direction ahead of the crack 
tip, and f(0) are known dimensionless angular functions (RICE, 1968). 

Under extensive creep conditions where the creep behavior of a material is described 
by the following steady-state creep equation 


(2) 


the asymptotic stress fields are given by 
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C* I/(n+ 1) 
6;,(r, 0) = (E : 6; (0). (3) 


Here, B is a temperature dependent coefficient in Norton’s creep equation, é,; the 
strain rate tensor, ¢, the von Mises equivalent stress and s,,; the stress deviator, while 
I, and 6,,(0) are dimensionless factors which are functions of n and 0, respectively, 
and tabulated by HRR (HUTCHINSON, 1968). 


2.2. Cavity nucleation and growth 


The cavity nucleation process is not yet fully understood, and here an empirical 
relation between the areal cavity density and creep strain is used to deduce the cavity 
nucleation rate as follows ; 


C = Aé, = F(T)o". (4) 


Here, C is the number of cavities per unit grain boundary area, C its time derivative 
and F(T) is a constant which is linearly proportional to the creep coefficient B. The 
relation has been used by several workers in creep rupture analysis and A is shown 
to be relatively independent of the test temperature (CANE, 1979; NEEDHAM and 
GLADMAN, 1980; RIEDEL, 1985). 

Compared to the cavity nucleation process, mechanisms of the cavity growth are 
better understood and cavity growth rates are derived when cavities grow by power- 
law creep, or by diffusion in a quasi-equilibrium or crack-like mode (HULL and 
RIMMER, 1959; CHUANG and Rice, 1973; RAJ and AsHBy, 1975; NEEDLEMAN and 


Rick, 1980). In the present analysis, we use a growth law proposed by TRINKAUS and 
ULLMAIER (1979) for the cavity growth by diffusion in a quasi-equilibrium mode. The 
volumetric growth rate of a cavity, V, is given by 


V=2nDo,,, where D= = (5) 
kT 
Here, o>, is tensile stress acting normal to the grain boundary, Dg is grain boundary 
diffusivity, d, is grain boundary thickness, Q is atomic volume, k is Boltzmann’s 
constant and JT denotes temperature. This law differs from others in that the cavity 
growth rate is not a function of the inter-cavity spacing, which is valid when the 
physical dimensions of elastic accommodations of grain boundary cavitation are of the 
order of the cavity size. When cavities grow according to (5), interactions among 
cavities with different nucleation times can be effectively neglected due to the non- 
overlapping internal stress fields built up by inhomogeneous deposition of atoms. 
The volumetric growth rate can be related to the areal growth rate of a cavity, A, 
by 


2 Fy 


=3 RY, (6) 


A 


where Fy and F, are geometrical factors describing the cavity volume and the grain 
boundary area occupied by a cavity, respectively. (Note that V = p*Fy and A = p’F, 
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where p is the cavity radius.) From (5) and (6), A under varying stress condition is 
written as 


ue 2 F 
A=k, | fos ar| >>, Where k, = “_ (2nD). (7) 
3 Fy’ 


For a cavity growth by creep, BUDIANSKY, HUTCHINSON and SLUTSKY (1982) 
showed that 


A = 2xBAo", (8) 
where x is a dimensionless factor given by 


~- 05 CF eae F069) 
AW) 


2no, n° 


and a, denotes the sum of the three diagonal stress tensor components. 


3. CREEP CRACK GROWTH MODEL 


3.1. A general solution 


Here, we follow the approach laid out by RrepEL (1983), and consider a stationary 
crack moving on the 0 = 0 plane after a certain incubation time (t,). Before the initial 
crack advance, a critical damage (a%,) is assumed to accumulate at a characteristic 
distance (Y,) from the crack tip, where X, is typically of the order of a cavity spacing 
at rupture. Under the continual cavity nucleation, the incubation time can be found 
from (TRINKAUS and ULLMAIER, 1979) 


(*t 


Lp aX. 0’) dt’ 


JO 


(*t 





ar| C(X., WH) A(X. t,t") dt, (10) 
0 


J0 


where C(r,t’), A(r,t’,t’), and a(r,t”) refer to the cavity nucleation rate at the time 
t’, the areal growth rate of those cavities at ¢” and the rate of increase of cavitated 
area fraction at f”, respectively, at a position r from the crack tip. 

After the incubation time a crack starts to move, and the damage at r = a—d)+ X, 
is the sum of the damage accumulated during the incubation period («,) and that 
during the subsequent crack growth («,), which are given respectively by 


a(a—ay+X,) = al C(a—a,+X., )A(a—a,+X., t,t”) dt’ (11) 
0 0 
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Fic. 1. Schematic diagrams of steady-state crack tip stresses and damage profiles. 


, “ da** ae * *\ 4 * * 7k* da* 
a,(a—ayot+ X,) = ar) J, C(a—a*+X., t®*)A(a—a**+X,, t*,t re * (12) 


Here, dy, a*, a** and a denote successive crack tip positions at initial time, /*, 7** 
and current time, respectively. Then, the equation of motion becomes 


Aap = a+. (13) 


The solution to this integral equation is a function of the applied load and the crack 
length increment (Aad = a—ay). After a crack moves more than r,,— X,, «; can be set 


to be zero in (13). Here, r,, is the nucleation distance on the grain boundary where 
local tensile stress reaches a critical value of the cavity nucleation (¢,,). 


3.2. A steady-state solution 


For a crack moving along the grain boundary (@ = 0 plane) with a constant velocity 
a, near-tip asymptotic fields and the profile of accumulated damage are invariant with 
respect to the crack-tip position as shown schematically in Fig. |. Let us take an 
arbitrary crack-tip position as a reference position and the corresponding time as a 
reference time (i.e. ¢ = 0), and assign r as a distance from the reference crack tip. 
Then, the amount of accumulated damage on the grain boundary at r = x (for x < r,) 
after a lapse of time ¢ since the reference time is given by 


t 


a(x, t) = a(x, 0)+| a(x,t,) dt, (14) 


( 


where the first and second terms on the RHS refer to a damage existing at r = x at 
the referred moment, and an additional damage accrued during the time /¢, respectively. 
Note that the first term in fact corresponds to the steady-state damage profile «(x). 
Under continual cavity nucleation, the damage accumulation rate at r= x at the 
moment ¢ = /, 1s given by 
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a(x, t;) = | | C(x, t) A(x, t, t;) dt, (15) 
0 


where C(x,t) is the cavity nucleation rate at the time t and A(x,t,1,) is the areal 
growth rate of these cavities at the time ¢,. For r > r,,, «(x,0) is zero because the local 
stress does not reach the nucleation stress. 


3.2.1. Crack growth under HRR fields 

3.2.1.1. Cavity growth by grain boundary diffusion. When crack-tip fields are HRR 
fields characterized by a steady-state value of the C* integral and cavities grow 
by diffusion, C and A in (15) can be respectively written, from (3), (4) and (7), as 


* |m/(n+ 1) 
C(x, t) = F(T)[é,(0)]” PS (x—at)~™/@r (16) 


A(x,1,t)) =k, | ‘(6 48 | ‘o(t)) 


C* I/(n+ 1) ]2/3 1/3 
40 (57) | (ats) sill 


*(x—dt)) Met Diy — at) @*  —(x—at,)"* "7 3 (17) 


Since the damage at x reaches «, when the crack tip moves to x— X,, a boundary 
condition to be used with (14) is that 


a(x,t) = %, when —_——, (18) 


Then, from (14), (16), (17) and (18), the steady-state damage profile can be written 
as 
a(x, 0) = a(x) = & —g(x, 4), (19) 
where g(x, a) is defined by 
g(x, d) = k,F(T)(a)~ PP C8Ct COD F(x) (20) 
with 


3 


1/2 
k = ke} | | [4.(0)]"(BI,)— 2+ mG 


f(x) oo : sis DL ame D(u" @+) _ yalet I)? 3 du. 


xX, 


An additional boundary condition to be used with (19) is 
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a(x) =0 when x=r,. 
Note that for x > X, and 5n—3m-+3 > 0, f(x) asymptotically becomes 


f( x) = Bim n)x' 5n— 3m+ 3)/(3(n+ 1)) 


where 


9(n+1)? 


; 25 
2n(5n—3m-+ 3) a 


B(m,n) = 


Subsequently, from (19)—(25), the steady-state crack growth rate d is found to be 
ele 
a _ [F( T)]° 5 fe Cre 3m)/(5S(n+ Tf (r,)]° 5 (26) 
R 


Here, r,, is obtained from (3) as 


where 


For r, > X., (26) can be rearranged as 


a =k,[F(T)]*°C*, 


where 


3/5 


k, |" 
k on fe | [B(m, n)]? *(k 5)" 3m + 3)/(S(n+ 1) (30) 


R 


Thus, a is linearly proportional to C* for any values of n and m which are constants 
describing the creep and cavitation behaviors, respectively. It is interesting to observe 
here that d@ is independent of X, values. 


3.2.1.2. Cavity growth by power-law creep. When crack tip stresses are still char- 
acterized by HRR fields but cavities grow by power-law creep, the areal growth rate 
of a cavity is given by (8) from which the following relation ensues ; 


A(x,t, t) = A, exp [22 | [o.(t)}" ai| 


* 


‘vain . 1/(n+ 1) + 4\ |/(n+ 1) 
BI (n+ 1) [(x—at) —(x—dt) ]/. 


(31) 


= A, exp |2«Bt6. (0)? ( 


Here, A, is the grain boundary area occupied by a cavity at cavity nucleation (i.e. 
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t = t). When the cavity nucleation stress exponent m equals the creep exponent n, the 
steady-state crack growth rate is approximately given by, from (15), (16), (31) and 
boundary conditions stated in (18) and (23), 


2x B[é.(0)}"(BI,) nlin+ Dn +1) 
a= 


For most cases, k,;C* » X, and a is found to vary linearly with C* as in the case of 
the diffusional cavity growth under C* fields. 


(C*)" O+D1(k,Ct)/ero_-y¥! vit: (32) 


3.2.2. Crack growth under elastic fields 

3.2.2.1. Cavity growth by diffusion.When a crack moves under the K, fields with 
a steady-state velocity ad, analyses of the creep crack growth can be carried out 
straightforwardly. Here, elastic stress fields given by (1) are used instead of the HRR 
fields, and respective substitutions of 1, 27, F(0) and Kj for n, BI,, &;(0) and C* of 
(3) can readily provide the corresponding relations. 

The steady-state damage profile equivalent to (19) is given by 


a(x) = op —ksF(T)(a)~ °° Ki?* w(x), 


ks= k i[f22(0)] > 51 F.(0))" (1/22)? * samp 


w(x) = \ 3u—”™ 2(/u—./X,)? * du. 
x, 


In (33), w(x) can be asymptotically expressed for x > X,. as 


w(x) = B(m)x®- >"> for m < 8/3, 


B(m) X87” for m > 8/3. 
Here, /(m) is a constant defined by 


18 
8—3m 


o(3)r(m-5) 


l(m—1) 


for m < 8/3, 


for m > 8/3, (39) 


where I is the Gamma function. 
As in the previous section, the steady-state crack growth rate can be deduced from 
(23) as 
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(KEEPS bP. 


Here, r,, is the nucleation distance given by 


r, =k.(K,)° where k, = BJ (41) 


2n| G, 


Using asymptotic expressions of w(x) given in (36)—(39), a can be rewritten as 
a=k,[F(T)]*°°(K)’ for m < 8/3, (42) 
k [F(T)}*°(K,) Cr? for m > 8/3, (43) 


where constants k; and ky are respectively given by 


ks + 3/5 3 
k, _ E [B(m)]’ Tk.) 3m)/10 


R 


3/5 


Ks 3/5 8 — 3m)/10 
kx = As | [B(m)]} °° [XJ OP". 


R 


Note here that ad ~ Kj for m < 8/3 and a@ ~ (K,)°°* > for m > 8/3. 


3.2.2.2. Cavity growth by power-law creep. When cavities grow by creep flow rather 
than by diffusion under elastic stress fields characterized by K;, the steady-state crack 
growth rate can be found analogously. For cases when the cavity nucleation stress 
exponent m equals the creep stress exponent n, a is found to be 


4x BLf.(0)]" 


2k Bap nl2 
In (: + F( me. (27)”"'~(n—2) 


a= Ki[x’? n)/2 — n) 3 (46) 


Note that ad ~ Kj for n > 2 and that ad ~ K? for n < 2. 


4. DISCUSSIONS 


4.1. Comparisons with experimental data 
4.1.1. 2.9Ni-1.1Cr steel. By conducting uniaxial creep tests with a 2.9Ni—1.1Cr 
steel, CHO (1984) obtained the following creep and cavitation behaviors at 823 K; 


é(s~') = 2.15 x 10> **[o,.(MPa)]’ (47) 
C(m~ 7) = 3.5x 10!%e. (48) 


These can be combined to give the steady-state nucleation rate of grain boundary 
cavities as 
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Fic. 2. Correlations between d and C* for a Ni-Cr steel studied by HONG (1989). Here, solid lines A and 

A’ are predicted @ values based on the cavity growth by diffusion at ¢, = 210° E and 5x 10°? E, 

respectively, while dotted lines B and B’ are predicted d values based on the cavity growth by creep at 
o, = 1x10-* Eand 5x 10° E, respectively. 


C(m~? s~') = 7.5x 10> ''[o,.(MPa)]’. (49) 


Creep crack growth rates of the same material are measured at the same temperature 
using the single specimen method and experimental details are described in Yu and 
HONG (1989). Experimental data presented in Fig. 2 are obtained from several 
specimens with varying initial loading conditions, and non-steady-state crack propa- 
gations are assumed to occur in the early transient response of the material. During 
that period the creep crack growth rate depends on the crack length increment (Aa) 
and applied load, but what we observe here is more or less a universal correlation 
between ad and C* which is typically expected out of the steady-state solution. Thus, 
a is reasonably characterized by C* including the non-steady-state crack propagation 
in this type of material, and this provides the basis of applying the present steady- 
state solution to the Ni—Cr steel data (HONG and Yu, 1989). 

Solid and dotted lines in Fig. 2 represent predicted d values from the models of 
the cavity growth by diffusion and creep, respectively, at two different levels of the 
cavity nucleation stress. The nucleation stress which best fits experimental data is 
2x 10° E for the diffusional growth and | x 10-* E for the creep induced growth, 
where E is Young’s modulus. According to the heterogeneous nucleation theory 
applied to the cavity nucleation (RAJ and AsHBy, 1975; Raj, 1978), o, can be less 
than 10-* E when the cavity geometrical factor F, is smaller than 10°. At present, 
there are no means to discern which nucleation stress is more realistic for ferritic 
steels doped with carbide particles. 

In the present analysis, the cavity growth equation proposed by TRINKAUS and 
ULLMAIER (1979) is utilized, which is strictly true when cavities are small and well 
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Fic. 3. Data by RieEDEL and WAGNER (1984) at 923 K are extrapolated to 1023 K using the reported 
activation energy. Here, solid and dotted lines are predicted ad values at 1023 K from the present models 
based on the cavity growth by diffusion and creep, respectively. 











separated. When inhomogeneities from neighboring cavities begin to overlap, depo- 
sition of atoms in the grain boundary becomes more uniform as in the rigid grain 
separation models (RAJ and AsHBy, 1975). In such a case, it can be seen from (29), 
(30), (21) and (7) that a is increased by (b/p)”° times without affecting the linear 
dependence of ad on C*, where denotes the inter-cavity half spacing. In a similar 
context, when the inter-cavity spacing is assumed to decrease uniformly with the 
continual nucleation of cavities as in ARGON (1982), the damage accumulation rate 
to be used in (14) is approximately given by 


10D 
h() 


which after all leads to another linear relationship between da and C* (YU and HONG, 
1989). Therefore, d increases in linear proportion to C* within the framework of the 
stress controlled cavity nucleation and the steady-state crack growth. 


a(x, t;) = [xC(x, t,)]}*?o(t,), (50) 


4.1.2. Nimonic 80A. From the experimental data by Dyson and others (1972, 1976) 
and Lovepay and Dyson (1979), the steady-state nucleation rate of Nimonic 80A is 
deduced as (FULLMAN, 1953) 


C(m-? s~') = 1.24x« 10> '®[o,,(MPa)]*. (51) 


In Fig. 3, experimental data by RIEDEL and WAGNER (1984) at 923 K are shown as 
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FiG. 4. Active cavitation zone developed ahead of the crack tip under the K, fields. 


solid circles, and those extrapolated to 1023 K using the crack growth activation 
energy of 64.4 kcal/mol (RIEDEL and WAGNER, 1984), are shown as open circles. In 
the calculation of d@ using (43), X, is taken as 6 y which is the inter-cavity spacing at 
rupture reported by DysON and MCLEAN (1972). Predicted a values based on the 
diffusional cavity growth show better agreement with the experimental data than 
those based on the cavity growth by creep. 

In (40), we have shown that a ~ K}***”°- [w(r,)]}°° and that w(r,,) is independent 
of r, for m > 8/3, which result in the insensitivity of ad on r,. Physically, this originates 
from the restrictive size of the active crack-tip damage zone which is much smaller 
than the nucleation distance as shown in Fig. 4. Since the active cavitation zone is 
virtually limited to several X,, r,, does not influence ad, and this seems to explain the 
relative difficulty of observing crack-tip cavities in superalloys compared to ferritic 
steels where the active damage zone usually spans more than 100 X.. 

For many materials which are known to have creep crack growth under elastic 
fields, experiments show that ad ~ K? where p is related to the exponent of the creep 
equation n. In section 3.2.2, we have shown for m = n that p = (2+3n)/5 (forn > 8/3) 
for the diffusional cavity growth and that p = n (for n > 2) for the cavity growth by 
creep. Experimentally observed p values are plotted against corresponding n values 
for various superalloys and a ceramic material in Fig. 5, and the model based on the 
diffusional cavity growth describes experimental data reasonably well. 


4.2. Comparisons with existing models 


In order to compare the present analysis with existing models in which aspects of 
the continual cavity nucleation are ignored, dependences of ad on C* or K, are presented 
in Table 1. For various situations of cavity nucleation and cavity growth, the following 
statements can be made: 
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Fic. 5. For various superalloys and a ceramic, experimentally observed p values (4 ~ Af) are correlated 

with corresponding n values (é, = Bo"). Respective data are for In-738 (NAzMy and WuUTHRICH, 1983), 

Nimonic 115, Astroloy and Inconel X-75 (FLOREEN, 1975, 1983), Nimonic PE16, Udimet 700 and Inconel 

718 (SADANANDA and SHAHINIAN, 1977, 1978), Nimonic 80A (RIEDEL and WAGNER, 1984) and for Alumina 
(BLUMENTAL and Evans, 1984). 


(1) Cavity growth by diffusion and crack growth under C* fields; The stress 
controlled nucleation model by Riedel predicts ad ~ (C*)'“"* while others 
show d ~ C*. 

(2) Cavity growth by creep and crack growth under C* fields; The critical strain 
model by Riedel predicts ad ~ (C*)""* ', while others show 4d ~ C*. 

(3) Cavity growth by diffusion and crack growth under K, fields; The model by 
WILKINSON and ViTEK (1982) predicts d ~ K; while the present model shows 


TABLE |. Dependence of a on load parameters for various models 





Crack Cavity WILKINSON 
tip growth and 
fields mechanism  VITEK (1982) RIEDEL (1983) Present model 


gZ a~ C* a~ ce DAqi” 1) a~ C* 


C* diffusion 


creep a ois C* d tee cr (n+ "(Aa' (n+ 0—¢,)t a oa c* 





g.b. . | a~ Kj?*?”?, m> 8/3§ 
diffusion a~ K7,m < 8/3 


a~ Kj,n> 


a~ Kj,n>2 
a~ Ky,n< 


' a~ Kj,n> 3t 
a~ K7,n<2 , ° 


creep 





* Critical strain model (RIEDEL, 1983). 
t Result by Hui and RiepeL (1981). 
§m: Stress exponent of cavity nucleation rate (~7). 
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TABLE 2. Activation energy of creep crack growth rate, Q,, for various 





models 
Crack Cavity WILKINSON 
tip growth and RIEDEL Present 
fields mechanisms Vitek (1982) (1983) model 
g.b. ) 7 
ne diffusion Cn @, Qe— n+l g. 2(Qa—@.)/5 
| 
ree 0 
creep n + ] 2. - ae AL, 3 
g.b. 
diffusion Or (208+ 3Q.)/5 
K, 
creep 0. OQ. d. 
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a ~ K\?**”> for m > 8/3. The latter seem to be the cases of experimental data 
shown in Fig. 5. 

(4) Cavity growth by creep and crack growth under K;, fields; All the models 
predict d ~ Kj forn > 2. 


One way to investigate the appropriateness of various models is to examine the 
predicted activation energy of creep crack growth rate, Q,, which is defined by the 


relation 
ad = a)exp  - | (52) 


Here, dy is a normalizing constant, and inherent Q, values out of various occasions 
are listed in Table 2. In deriving Q,, note that B and D are the only temperature- 
dependent material parameters which contain activation energies of creep (Q.) and 
grain boundary diffusion (Q,), respectively. For most ferritic and stainless steels 
which are tested at 0.4—0.5 7, where the activation energy for creep is close to that 
of grain boundary diffusion (ROBINSON and SHERBY, 1969 ; SHERBY and YOUNG, 1975; 
Nam, HonG and SHIN, 1983), Q, is found to be very small which is consistent with 
most models except the one by RiepEL (1983) which treats the diffusional cavity 
growth under C* fields. Similarly, for the crack growth under K;, fields, present model 
based on the diffusive cavity growth predicts the Q, value of a Nimonic 80A reasonably 
well because values of Q,, QO, and Q,. of the same material are reported to be 64.4, 
27.4 and 72-96 kcal/mol, respectively (ASHBY, 1972; SMITHELLS, 1983; RIEDEL and 
WAGNER, 1984). When cavities are assumed to grow by creep, all models show that 
QO, = Q., which tentatively agrees with experimental results. 


5. CONCLUSIONS 


For macroscopic crack advance by the nucleation, growth and interlinkage of 
grain boundary cavities, steady-state crack growth rates are derived based on the 
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assumption that the cavity nucleation occurs at the attainment of a critical stress (¢,,) 
and continues in a strain-controlled manner. Main results are summarized as follows; 


(1) When the crack growth occurs under HRR fields, ad depends on oa, and varies 
linearly with C* regardless of the cavity growth mode. Analyses of a Ni-—Cr steel 
suggest that a, lies between 10° *-2 x 10° E. 

(2) When a growing crack is embedded in elastic fields, d ~ K{°**”’° (n > 8/3) for 
a cavity growth by grain boundary diffusion and ad ~ K{ (n > 2) for a cavity growth 
by creep. For both cases, d is independent of o,, and depends on X, which is about 
the order of intercavity spacing at rupture. A collection of data on superalloys can be 


better described by the former. 

(3) Size of the active cavitation zone ahead of a moving crack tip is much smaller 
under elastic fields than under HRR fields. 

(4) Activation energies of creep crack growth predicted by the present model tend 
to agree reasonably well with those of a superalloy and various steels. 
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ABSTRACT 


THE DIFFUSIVITY of ultrasound in an untextured aggregate of cubic crystallites is studied theoretically with 
a view towards nondestructive characterization of microstructures. Multiple scattering formalisms for the 
mean Green’s dyadic and for the covariance of the Green’s dyadic (and therefore for the energy density) 
based upon the method of smoothing are presented. The first-order smoothing approximation used is 
accurate to leading order in the anisotropy of the constituent crystallites. A further, Born, approximation 
is invoked which limits the validity of the calculation to frequencies below the geometrical optics regime. 
Known result for the mean field attenuations are recovered. The covariance is found to obey an equation 
of radiative transfer for which a diffusion limit is taken. The resulting diffusivity is found to vary inversely 
with the fourth power of frequency in the Rayleigh, long wavelength, regime and inversely with the 
logarithm of frequency on the short wavelength, stochastic, asymptote. The results are found to fit the 
experimental data. 


NOMENCLATURE 


A, B,C Stokes-like parameters for the radiation intensity, (6.7) 
Ao, Bo, Ai, By expansion coefficients for A and B, (6.13) and (6.14) 
rotation matrix element, (2.7) 

stiffness tensor, (2.1) 
ensemble averaged stiffness, also the Voigt estimate of effective stiffness, 
(2.2) 
effective (Voigt average) longitudinal and transverse wave speeds 
generic wavespeed 
diffusivity, (1.1) 
diffuse field energy density in transform domain, (4.5) 
unaveraged Green’s function, (2.1) and (3.2) 
Green’s function of the bare, Voigt average, medium, (3.3) and (3.7) 
ensemble average Green’s function, (3.3) 
bare longitudinal and transverse propagators, (3.10) 
mean field propagators, (3.22) and (3.23) 
covariance of the Green’s function, (5.1) 
Latin Cartesian indices 
unit dyadic 
diffuse field energy flux in transform domain, (4.16) 
ratio c,/c; of effective medium wave speeds, (7.6) 


+The author carried out this work while on leave at the Department of Mechanical Engineering, 
University of California, Berkeley, CA 94720, U.S.A. 
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k intensity operator, (5.4) 
L(0), M(0), N(@) Eq. (2.15) 
self energy operator, (3.3) 
wave vectors 
unit wave vectors 
mean field energy propagators, (6.1) 
covariance source, (5.7) 
dimensionless measures of frequency, (7.4) 
configuration space position vectors 
generic attenuation 
longitudinal and transverse wave attenuations, (3.29), (3.30) and (1.1) 
Opp, Xer, Xrp, Arr partial attentuations, (3.35) 
A ee oe angle weighted partial attenuations, (3.35) 
Greek Cartesian indices 
microscale inverse length, (2.5) 
fluctuating modulus, (2.2) 
double mean field propagator, (5.5) 
Spatial Fourier transform parameter for mean square quantities, (4.5) 
Kronecker delta 
three-dimensional Dirac delta function 
infinitesimal positive constant, (3.2) 
microstructural two point correlation function, (2.4) 
Fourier transform of n(r), (3.16) 
Eq. (3.33) 
the third independent modulus of a cubic crystal, (2.6) 
expansion coefficients for R’ and R’, (6.15) 
self energies for the mean field propagators, (3.14) and (3.19) 
generic self energy 
wave frequency 
temporal Fourier transform parameter for mean square quantities, (4.2) 
constant eighth rank tensor, the dimensionless modulus variance, (2.4) 
0/0,, the partial spatial derivative 
complex conjugate, (4.3), also @ ~ w+Q, (5.1) 
ensemble average 


1. INTRODUCTION 


THE MICROSTRUCTURE Of polycrystalline aggregates has been probed by means of 
ultrasonic scattering for at least 40 years, since the work of MASON and McCSKIMM 
(1947). That ultrasound with frequencies of the order of 10 MHz and inverse wave 
numbers of order 50 microns would be sensitive to microstructural variations on the 
length scales of relevance in typical polycrystals has long been appreciated. The names 
of PAPADAKIS (1968, 1981), BHATIA (1959), Adler (FittING and ADLER, 1981), 
Goebbels (GOEBBELS and HOLLER, 1978 ; GOEBBELS, 1980; Guo et al., 1985) and VARY 
(1978) are amongst those associated with laboratory efforts in the pursuit of the 
ability to assess microstructure, and potentially therefore, mechanical properties, using 
nondestructive ultrasound. In spite, though, of considerable effort over many years, 
the most common technique, that of ultrasonic attenuation, though proven in the 
laboratory, has not seen much application in the field. 
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That technique calls for the measurement of the exponential rate of spatial decay 
of an ensemble averaged plane wave. Due to a host of potential systematic errors, the 
unambiguous measurement of that rate is, however, quite difficult (TRUELL et ai., 
1969). The finite beam width to wavelength ratio of ultrasonic transducers is responsible 
for geometric, non exponential, decays as the beam transmits from near to far field 
conditions. Most studies are carried out in slab-like geometries and the wave amplitude 
monitored as it reflects between opposite faces. Imperfections in the reflection 
coefficients at these faces can artificially enhance the apparent attenuation. Faces 
which are insufficiently parallel can distort the apparent pulse amplitudes. Internal 
friction, or absorption, which is more properly viewed as a temporal decay than a 
spatial one, contributes in such a way as to be indistinguishable, in these configur- 
ations, from scattering based decay. High attenuation rates are immeasurable except 
over distances so short as to bring the ergodic hypothesis, crucial to comparisons 


\ 


between theory and experiment, into severe question. m 
In part due to the above difficulties there has been an increasing effort in recent 


years to measure microstructure from the energy which is incoherently and singly 
backscattered from a beam (GOEBBELS and HOLLER, 1978; Fay, 1973; SANIIE and 
BILGUTAY, 1986; SANtIE et al., 1988). The potential for systematic error in these 
configurations is perhaps less severe than in the conventional parallel wall slab con- 
figurations. The technique is also attractive in that it makes fewer demands on 
specimen geometry and can potentially give information on microstructure as a 
function of depth. The chief difficulties appear to be traceable to the need for extensive 
spatial averaging of signal powers in order to minimize the effect of the inevitable 
meaningless fluctuations of an incoherent wave field, the need for minimum slab 
thicknesses in order to resolve backwall echoes and distinguish them from the import- 
ant incoherent grass and the absence of a complete theory relating backscattered 
power to microstructure. 

It is proposed here that measurement of the evolution of the incoherent and 
multiply scattered, fully diffuse, ultrasonic wave field in a material with a scattering 
microstructure may provide a more robust measure of that microstructure. The 
published literature addressed to such wave fields is meagre, and seems to be confined 
to the paper by Guo et al. (1985) and to two by WEAVER (1989a, b). WEAVER (1989a) 
considered sub MHz range ultrasound in a specimen with an artificial centimeter scale 
microstructure. While not directly applicable to the usual materials of interest in 
NDE, it did demonstrate, over several orders of magnitude, the applicability of the 
concept of diffusion to the evolution of these fields. That is, it was shown that the 
ultrasonic energy density evolved in accordance with a diffusion, or heat, equation, 
modified by the inclusion of an extra term representing temporal decay. 


0E/ét = DV’ E—€E+source (1.1) 


where E is the elastodynamic spectral energy density at frequency w; D(q) is the 
frequency dependent diffusivity and (q@) is the absorption rate, taken to be uninter- 
esting for the purposes of the present communication and henceforth neglected. 
WEAVER (1989b) has also shown strong correlations between measured diffusivities 
in steels and their heat treatments and fracture toughnesses. Guo et al. (1985) after 
a study of such fields in a range of materials, also report a good fit to the predictions 
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of a diffusion model. Guo et al. and Weaver have both emphasized that the diffuse- 
field technique promises to be able to distinguish between temporal decay, or absorp- 
tion, and scattering. Each is potentially a nondestructive testing parameter of interest, 
the former having received little attention in the past due to the difficulty of its in situ 
measurement, (e.g. TVERDOKHLEBOV, 1986). 

No theory for the diffusivity of an ultrasonic field in a polycrystal is immediately 
obvious. Guo et al. (1985) suggest, in analogy to a random-walk model, that the 
diffusivity D should be one-third of the product of a mean free ray path and the wave 
speed. They then identify that mean free path with the inverse of the attenuation and 
suggest that the diffusion constant D is simply an inverse measure of attenuation. The 
argument is attractive due to its simplicity, but it must be modified to account for the 
presence of different wave modes, the transverse (T) and longitudinal (L), each with 
its own wave speed c and attenuation «. We must further recognize that the appropriate 
mean free path should be half of the inverse of the attenuation (as it is the energy 
which is diffusing and the attenuation rate of the coherent energy is twice that, «, of 
coherent amplitude). With these identifications we write, as a conjecture, 


D = fr(Cr/6ar) + fi (C,/6%,), (1.2) 


where f; and ff, = 1—/; are the mixing weights for the transverse and longitudinal 
contributions to the diffusivity. Their values are not obvious. The only a priori likely 
possibility is 


Self = 2ci/c7 (1.3) 


corresponding to the equilibrium partition of the energy of a fully developed diffuse 


wave field (WEAVER, 1982). 

Equation (1.2) has the virtue of being dependent only on quantities available from 
standard theory and/or experiment. As Guo et al. indicate it also implies that the 
diffusivity should, for frequencies below the geometric optics regime, scale inversely 
with frequency with a power between two and four. This follows from the standard 
theory (for which the reader is directed to the recent work by STANKE and KINO, 1984 
or HIRSEKORN, 1982, 1983, 1985, 1986 or KARAL and KELLER, 1964) or to the general 
discussions (PAPADAKIS, 1981; BHATIA, 1959; FiTTING and ADLER, 1981) for the 
attenuation of mean waves in polycrystals. That theory indicates that attenuation 
scales with the fourth power of frequency in a domain for which wavelength is much 
larger than microscale length (the ‘Rayleigh’ domain) and with the square of frequency 
in the opposite limit (the ‘stochastic asymptote’) but is ultimately bounded by the 
geomtrical optics limit where «~' is of the order of the microscale length. 

The frequency dependence of the diffusivity observed by Guo et a/. was however 
very clearly an inverse first power. In consequence, (1.2) cannot be correct. We will 
see in this communication that that error is due to the identification of that “mean 
free ray path” required by the random walk model with 1/2«. The model requires a 
mean free path which corresponds to the distance travelled by a typical ray before it 
is significantly scattered away from its original direction. At low frequencies, in the 
Rayleigh regime, grain scattering is isotropic; equal energies are scattered in all 
directions ; hence 1/2« does correspond to a mean free path. At higher frequencies 
the scattering is increasingly biased towards the forward direction; typical rays are 
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scattered through small angles; hence the appropriate mean free path is somewhat 
greater than 1/2a. 

The present communication is addressed towards the derivation of a more accurate 
expression for D. The derivation is confined, like the attenuation calculation of STANKE 
and KINO (1984), to the case of an untextured aggregate of cubic-symmetry crystallites. 
This is the simplest polycrystal and includes iron as a special case, allowing comparison 
with the results from one of the specimens studied by Guo et a/. After discussion of 
the assumed statistics of the polycrystal in Section 2, Section 3 presents, by means of 
consideration of the ensemble average Green’s function the standard result for the 
attenuations, explicit expressions for which are given in Section 7. Sections 4 and 5 
present the analogous theory for the mean square Green’s functions, and therefore 
for the energy. Section 6 solves the resulting equations of radiative transfer and 
concludes with an expression (6.26) for the diffusivity. That diffusivity is found to be 
given by a function of longitudinal and transverse wavespeeds and attenuations, as 
anticipated in (1.2), but to depend upon scattering-angle weighted versions of the 
attenuations as well, as anticipated above. Upon choosing a form for the micro- 
structural spatial correlation function, the diffusivity is evaluated in Section 7. It is 
shown that the frequency dependence of D spans a range from strong to weak. The 
low frequency dependence is with the inverse fourth power of frequency. At high 
frequencies but below the geometric optics limit D~ ' scales only logarithmically with 
frequency. Section 7 concludes with recommendations for further work. 


2. MATHEMATICAL PRELIMINARIES 


Green’s function (dyadic) for the response of a heterogeneous anisotropic passive 
elastodynamic medium is given by the causal solution to 


[ — 507 + Oj Cyuij(x)0 JGin(x, x)= 6,0°(x—x’)d(2). (2.1) 


In (2.1) units have been used such that the material density is equal to one, and 
G;;(x, x’) is the displacement response in the i direction at position x due to a unit 
impulse applied in the / direction at position x’ at time zero. 

The usual fourth rank stiffness tensor is given by 


Crirj(X) = Cheng + Yeirj(%) (2.2) 


<y(x)> = 0. (2.3) 


Angular brackets represent ensemble averages, hence C has a mean value of C° and 
y represents the modulus fluctuations. We shall be considering the system (2.1) for 
the case y small, and discussing mean responses and mean square reponses to leading 
order in the magnitude of y. As such, all relevant statistical information regarding the 
material heterogeneity is contained in the covariance function 


Pair (X)Vapys(X’)> = Beir n(x —x’|). (2.4) 
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Any other statistical information about the medium is necessarily of order y” with 
n > 2, and therefore negligible. 

Two assumptions regarding the statistics of the polycrystal are implicit in the form 
of (2.4). The first assumption is that the tensorial character of the covariance function 
decouples from the spatial dependence. This is equivalent to an assertion that there 
are no orientation correlations between different crystallites. The assumption is a 
standard one. We have also assumed that the polycrystal’s second order statistics are 
homogeneous (a function of x—x’ rather than x and x’ separately) and isotropic 
(independent of the direction of x —x’). 

If v is taken to be a dimensional measure of y then we define n(r) as v times the 
probability that two points separated by a distance r lie within the same crystallite. 
Hence the eighth rank tensor = is dimensionless. Like STANKE and KINO (1984), we 
shall choose an y of the form 

n(r) = v?e”’. (2.5) 


The present work will not, however, call for any specification of 4 until necessary, in 
Section 7. We will occasionally cite B~', though, as a measure of microscale length. 
With respect to crystal axes, the moduli of a cubic crystallite are 
Cixi = M5041 + MOO jp + 500 jx) + VO; jx: 
= Cijxi - VO i jks (2.6) 
where the last term in each line is not a tensor, and is defined as vanishing unless all 


indices are equal, and being unity otherwise. 
With respect to laboratory axes the crystallite modulus tensor is given by 


3 
Cis = Chics tv py GQ} A. 


n= 1 


= Chins tijnt (2.7) 


where a? is an element of the rotation matrix representing the transformation between 


the crystallite axes and the laboratory axes. 
If the rotation between crystal and laboratory axes is represented by the three Euler 


angles @, © and ¢, the rotation matrix elements are given by 
a; = —cosOsindsinf+cosdcosl, a} =cosOcos dsin(+sindcosf, 
; = sin sin O, 
= —cosOsindcosf—cos@sinl, aj; =cosOcosdcos(—sing sin, 
cos ( sin ©, 
=singdsinO, a;= —cosdsinO, a} =cos®, (2.8) 


as found in Eq. (2.50) of Bunge’s treatise (BUNGE, 1982). 

We now further assume that all crystallite orientations occur in the ensemble with 
equal weight. This corresponds, when taking averages of quantities which depend on 
a, to the use of the “uniform measure” (BUNGE, 1982) 


sin@/8n7dOdodt [0<O<n, 0<o<2n, 0<6<2z] 


when integrating over crystallite orientations. 
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The ensemble average modulus of our aggregate then becomes 


: inOdOdd¢d 
(Cini? = Cis = Cy aapata = p ¢ 





= Cini + «tiji>: (2.9) 


The tensor <t)> is manifestly an isotropic fourth rank tensor which is invariant under 
any permutation of its indices. There is only one such tensor 


<tijni> - (5; bx: + 54.5 jp +56 %)|1), (2.10) 


where the magnitude |f]} may be determined by evaluating any one component of 
(2.9). For example <73333> is, according to (2.10), equal to 37). It is also, according 
to (2.9) and (2.8), equal to 3v/5. Hence |t| = v/5 and 


Ctijcr> = (¥/5)[6;Ox0 + 5.6 yp +50). (2.11) 
Therefore the Voigt-average, and isotropic, stiffness C° is given by 
Pict = (A! + ¥/5) 5; De: + (u! + v/5) [646 jp +56 x] 
= (€7 —2¢F)6 Ox: + C76 x6 jp + 5:6 x] (2.12) 


where c, and c; represent the longitudinal and transverse wavespeeds in the Voigt- 


average medium. 
The variance, (2.4), is an eighth rank tensor: 


va = (Ciupys Cijxi> ” (Cups >< Cijnr> 
= Cbapystijzi> — (tapys>< tiers (2.13) 


where each of the factors of C is evaluated at the same point, in accordance with (2.4) 
evaluated at x = x’. This tensor is manifestly invariant under permutation of Greek 
or of Latin indices. It is also invariant under exchange of all the Greek for all the 
Latin indices. It is isotropic and not a pseudo-tensor, it is therefore constructable 
entirely from Kronecker deltas. Correspondingly, each index of any nonvanishing 
component must occur an even number of times. There are three independent tensors 
with these properties. The general linear combination is 


<tt> —<t><t> = bv? {bi + 5x6 jp + 5116 jx} {528555 + 52595 + 5555, 

+ dv? {6,i5p 6,55, plus all permutations consisting 
of exactly two pairings between Greek 
and Latin indices, one pairing between 
Greek and Greek indices, and one pairing 
between Latin and Latin indices, 72 terms in all} 

+ hv? {6,;6;6,,65; plus all permutations consisting 
of exactly four pairings between Greek 
and Latin indices, 24 terms in all}. (2.14) 


The numerical coefficients b, d and h may be determined by evaluating three 
independent components of =. For example 
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= t3222l2333> — (t3222)<l2333> 
= 9dy’ 
= 6v*<(a3)(a3)* (a3)(a3)*> + 3v?<(a3)*(a3)*> 
—v?/70. 
Hence 
d= —1/630. 
The other coefficients may be found by similar evaluations. They are 
b = 13/315—1/25 
and 
h = 1/180. 


Certain inner products on & will be required for the present calculation. They are 
tabulated here. In terms of two unit vectors p and §, with angle 6,, between them, they 
are 


l 4 
595 0°S ? On. + 595 09S é.. 


M(6,;) = Sin! PaPiSpS,P, Pr, 


= 595 + 595008 Op 


N(O,,) = _ EM bi PpPjSy §,05; = 


Py a. 29... 
"os as" * 


3. ENSEMBLE AVERAGE RESPONSES 


The bulk of the theoretical work on waves in polycrystalline media has focussed 
on attempts to describe an ensemble average plane wave or more specifically, the 
effective wave number of the average field. This wavenumber is normally taken to 
indicate the speed and the attenuation of waves observed in the laboratory. The 
present section presents a derivation equivalent to the usual discussions, but in a 
language more suited to the requirements of later sections. It concludes with closed 
form expression for the average Green’s dyadic in the frequency domain and 
expressions for the attenuations which are equivalent to those derived elsewhere. The 
focus is on the mean response, which is calculated within an approximation, variously 
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called the first-order-smoothing approximation or Keller approximation, valid to 
lowest order in the magnitude, v, of the moduli fluctuations. The final results are 
expressed in closed form valid for the frequency domain below the geometrical optics 
limit. 

The starting point is (2.1) together with the (temporal) Fourier transform pair 


Goin = | G;,(t) dt, 


| x 
G;,(t) a a | Grin e e dw. (3.1) 


The causality condition or alternatively, the quiescent initial conditions, together with 

the boundedness of the energy, assures that G,, is analytic for all w in the upper half 

complex plane and in particular for all @ with infinitesimal positive imaginary part. 
In the w-domain the governing PDE for the Green’s dyadic becomes 


{(@+ ie)" 5, + CpijOn0 + On Yeu j(X)O j} G(X, x)= 5,,0° (x —x), 3.2) 


where an ie has been included to emphasize the infinitesimal positive imaginary part 
of w. The w subscript on G has been suppressed as well. It will be re-included at 
occasional future points in the present exposition, but only when necessary for clarity. 

Equation (3.2) is of a form considered by others (KARAL and KELLER, 1964; FRISCH, 
1968 ; McCoy, 1981). The last term in brackets is a stochastic operator. The assumed 
smallness of y suggests the use of an expansion in powers of y. But as FRISCH (1968) 
emphasizes, any simple perturbative expansion of G in powers of y will fail at large 
|\x —x’|. The difficulty is overcome in most expositions by the method of “‘smoothing” 
described in Frisch’s case very elegantly using diagram methods. The result is an 
integral equation for <G), the Dyson equation 


(GiAX,X’)> = G(x, x) + | | Gig (x, yp (¥,2)<G,(z,x)> By d?z (3.3) 


where m is the ““mass”’, or “self-energy” operator and G° is the known Green’s dyadic 
for the bare, Voigt-average, medium, i.e. the solution to (3.2) when y = 0. It is given 
below. 

The Dyson equation is exact. An approximation is, however, necessary for the 
evaluation of m. The “First-Order-Smoothing Approximation” (FOSA) (FRIscH, 
1968) approximates m by its lowest order term in an expansion in powers of the 
fluctuating moduli. It is generally argued that FOSA is exact in the limit of small y. 
This corresponds in the present case to sufficiently small values of a non-dimen- 
sionalized |v|.¢ The FOSA expression for m is 


+The important issue as to how small must |v| be for FOSA to be accurate will not be completely 
addressed here. In the case of iron, |v|, when measured in the most obvious nondimensional manner, |v|/ 
c7, has a value of 1.66. On the other hand, STANKE and KINO (1984) argue for a nondimensional |v| in the 
case of iron of order 0.14 or less. 
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0 0 0 0 
mg (Yy, Z) ~ (4 Yapys(Y) ay. Gi (y, Z) a, Vika(Z) = ) (3.4) 


Vs zZ l 


The Dyson equation (3.3) may be solved with ease in the domain of spatial Fourier 
transforms. We define the spatial Fourier transforms of G°, m and <G). 


l 


Gi,(p)d*(p—q) = ays | [a'xaxe ”*Gi,(x,x’)e*™, (3.5) 


where the diagonal, or delta-function, character of G° in transform space follows 
necessarily from the homogeneity of the bare medium. Equation (3.5) has an inverse 


l | is 
G%x) = Go | | d’pd°qe”*Gi,(p)5*(p—qye*™. (3.6) 


G°(p) is readily constructed by means of the Fourier transform of its governing 
equation, (3.2) with y set to zero. 
{(@ + i€)*51 — PP jCeiij} GzP) = Sia- (3.7) 
In direct notation the equation for G°(p) becomes 
(pp{ (w+ ie)’ —p*c7} + I—pp) {(@+ ie)’ —p*c7})°G° = 1. (3.8) 
The solution to (3.8) is available by inspection: 
G°(p) = ppg’ (p) + (I—pp)g’"(p), (3.9) 
where g*” and g’’ may be termed the bare longitudinal and transverse propagators 
9°" (p) = [(w+ie)?—p*ci]~', 9°" (p) = [(@+ie)* —p*c7]'. (3.10) 


It is of particular interest to note that, for infinitesimal positive ¢, the imaginary parts 
of these propagators are given by 


Img” (p) = —nsgn (w)d(w* —p* cj 


Img?’ (p) = —2sgn (w)d(w? —p*c7). (3.11) 


The Fourier transform of <G) is also diagonal by reason of the statistical homo- 
geneity of the medium. 


<Gi(P)>d°(P—Q) = on? | [exarw e®*<Gin(x, x’) e*™, (3.12) 


(Giz (x, x’)> ” 


| |ervarver<c.@nrer—ae ee (3.13) 


(2n)’ 


The Fourier transform of m is of the form 
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l — 
—o,,(p)d°(p—q) = On) {ors d?x’e” ®*m, (x, x’) e*™. (3.14) 


Like <G), m is diagonal in the Fourier domain in consequence of the statistical 
homogeneity of the medium. 
Substitution of the FOSA approximation for m (3.4) into (3.14) yields 


O5)(P) = PuPi fe rGx (8) Ein? 0,0;{n(r) e-”*}. (3.15) 


It is more convenient to have o expressed in terms of G° evaluated in the Fourier 
domain. Hence we define a spatial Fourier transform of the two point correlation 
function 


i(p) = [arr e-™ (3.16) 


l 
(2n)* 
and find that (3.15) may be re-expressed as a convolution between G° and 

05(P) = —P.P1 [a°scatsmie—s, sit" 
In direct notation 


a(p) = - [a's Go .E i (p—s). (3.18) 


By general symmetry considerations o must be, like G° and any second rank tensor 
function of a single vector p, of the form 


o(p) = a (p)pp+ oa’ (p)(1— pp). (3.19) 


Comparison of (3.18) and (3.19) with substitution of (3.9) gives expression for o” and 


o*(p) = — fares {39°"(s) + (I—3)9°" (s)}#(P—s) (3.20) 


| ; . : 
o'(P) = —5 | d*s(I— 8) pe. {39°"(s) + (I—5)g”’ (s)}i(Pp—S). (3.21) 


A spatial Fourier transform of the Dyson equation and an invocation of the 
definitions (3.5), (3.12) and (3.14) reduces the Dyson equation to an algebraic equation 
for the Fourier transform of <G)». The solution is 


<G(p)> = 9° (p) bP +9" (v-)A— pp) (3.22) 
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g' (p) = (@*? —p’ci+o'(p)) | (3.23) 


g' (p) = (w@? —p’cz +a" (p))'. (3.24) 


These are the propagators for the mean field. They reduce to the bare propagators 
g’’ and g’’ when vy, and therefore o, vanishes. The term in epsilon has been dropped 
as the small but nevertheless finite imaginary parts of the os will dominate the 
infinitesimal e. 

The wavenumbers of the mean field are extractable directly from (3.23) and (3.24) 
by seeking that p at which the g are singular. Were one to take the inverse Fourier 
transform (3.13) of the expression (3.22), the major contribution would arise from 
the pole near |p| = w/c. The effective wave numbers are therefore the solutions, p, to 
the algebraic conditions 


w*—p*ci +o*(p) =0 (3.25) 


w*—p*c7+a'(p) = 0. (3.26) 


Evaluation of the effective wavenumbers requires the analytic or numerical solutions 
to these equations where the os are given by (3.20) and (3.21). The solutions simplify, 
however, in what is sometimes called the Born approximation (STANKE and KINO, 
1984). That approximation is motivated by a recognition that the os scale with v* and 


are therefore small. The solutions, p, to 3.25 and 3.26 are therefore close to the bare 
solutions p = w/c. We are thus led to approximate o(p) with o(w/c) and to write 
closed expressions for the effective wavenumbers 


p'! = (w*/cz +0" (@/cz)/c7)"’, (3.27) 
pr! = (w’/ci +0" (w/c,)/cz)"”. (3.28) 


As KARAL and KELLER (1964) discuss, this approximation fails in the very high 
frequency, geometrical optics, regime. The approximation is by no means a necessary 
one and STANKE and KINO (1984) give numerical solutions to the more general 
equations (3.25) and (3.26). We, though, retain the approximation because it simplifies 
the later analysis and is in any case valid except at very high frequencies. 

The attenuations are the imaginary parts of the effective wavenumbers and are 
given, to the same approximation, by 


a? = Imo’ (w/c,7)/2@cr (3.29) 


a” = Imo*(w/c,)/2@c,. (3.30) 


The imaginary parts of the os are in turn determinable from (3.20) and (3.21) together 
with the identification (3.11). 





Diffusivity of ultrasound in polycrystals 


a 
Im a! (jw Cc; ) = > 


@ 


T ie w° oOo .@ 
— ps — —— ee a atl ail 2a 
+ ‘lL. (I—S)G(1—j):--S P i (a: - )a , (3.32) 


T T 


where the integrals over the magnitude of the vector s have been absorbed by the 
delta-functions in (3.11), leaving integrals over the unit sphere 8. 

In terms of the three functions of the angle @,, between the unit vectors p and § 
which are defined by 


0 (Ops) = H(po/c, —Se/c,), 
n’" (0,5) = n'*(0,.) = H(po/c, —S/cr), 
n'" (0,.) = i(po/er —$/cr), 


L I 
OSA, 4 = Apr t+ %z,, 


+ 


n?>w* 


2c} l 


4 
7 n'(8)L(8)d cos 8, 
L 


=— | din“ (0,,)L(0,,) = 


4 
. 


3 


| ds"? (0,,){M(O.,) — L(O,y)} 


4 l 
, | n'" (0)(M(0) — L(0)) dcos 8, 
303 |, 


9 
47, = 5} (Cr/C,.) 4x7, 


nmu* 


T 


| d?sn7" (0,,){N—-2M +L} 


4 +1 
= | n™™(0)(N(0) —2M(0) + L(0)) dcos 0. (3.35) 
T 1 


Four additional attenuation-like quantities may also be defined, «,, a7, a, and 
x“, by expressions identical to those of (3.35) except that one extra factor of cos @ is 
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inserted in each integrand. The primed quantities appear to represent the preference, 
if any, of the scattering for the forward direction. It is worth noting that in the Rayleigh 
limit w « fe, all the ys become independent of their arguments. In consequence, all 
attenuations scale with w*, and all primed attenuations vanish. 

Upon specification of n(r), the above integrals can be evaluated. For the case of an 
exponential n(r) (2.5) this is done in Section 7. 

The major qualitative difference between the bare medium and the effective medium 
is the presence of this attenuation in the effective medium. It should be clearly 
understood, though, that this attenuation is not due to a dissipative process. In a 
sense it is entirely an artefact of having asked for the ensemble average response. 
Phase differences accumulating between the disturbances in different realizations lead 
to destructive interferences upon averaging over the ensemble. Viewed alternatively, 
the attenuation described above can be thought of as due to the scattering of energy 
out of a beam. In any case, there is no true dissipation present. Therefore the energy 
which is associated with the attenuating mean field, itself attenuating also, is a very 
poor representation of the actual energy in a typical member of the ensemble. This 
latter quantity is probably more accurately described by the mean of the energy, not 
the energy of the mean. In order to describe experiments such as that of Guo et al. 
(1985) or WEAVER (1989a,b) or backscatter studies (GOEBBELS and HOLLER, 1978; 
Fay, 1973; SANE and BILGUTAY, 1986; SANIIE et al., 1988) all of which scrutinize 
mean square responses, a different theory is required. The next section defines a few 
energy-related quantities of interest and the following section presents the multiple- 
scattering theory for their evaluation. 


4. ELASTODYNAMIC ENERGY, FLUX AND DIFFUSIVITY 


The energy density at position x and time ¢ due to a unit impulsive source at x’ 
polarized in the j direction, and for simplicity summed over that direction, is given 
by the sum of kinetic and strain energy terms 


e = 1/2G;,,(x, x’, )G,,(x, x’, 2) + strain energy density. (4.1) 


By means of the virial theorem we may recognize that the total energy density is equal, 
on suitable temporal or spatial averaging, to twice the kinetic energy density. Thus 
we are led to define a quantity which appears to be a temporal Fourier transform of 
the energy associated with a band passed process, 


éQ,x,x)=Y] ec ™|BR)G,(x, x’, |? de, (4.2) 
ij J—-@® 
where B represents a moderately narrow and moderately well damped band-pass 
process, centered, typically in polycrystalline applications at a frequency of order 5 
to 50 MHz, at central frequency w with width Aw. Q is of the order of the inverse of 
the time scale upon which the relatively slow evolution of the diffuse energy is to be 
studied ; Q « Aw « w;Q « ac « w and (X) represents a temporal convolution. 
By employing the Fourier representation (3.1) for G(t), é may be re-expressed as 
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dida dev’ e~ au - it + ie’ t , 
é= we Qn)? B(w) B* (w’ Jaw’ Gij(X, x ‘Gt ij(X, x ), 


(4.3) 


where an asterisk denotes a complex conjugate. The time integration yields a delta 
function which then allows ready evaluation of the w’ integration, and é becomes, in 
direct notation, 


é= \e y Bo) B*(o+ Qao(w + Q)G, (x, x’): G*, o(x, x’). (4.4) 


As B(q) is a smooth function with width Aw > Q, and as w » Q, 


é dw mao , 
(Qn 7a) | BCa)|?e w°G,,(x, x’): G3, o(x, x’). 


The above equation indicates that é is a frequency weighted average of the product of 
Green’s function at two different frequencies. By an ergodic hypothesis this frequency 
average is equivalent to an ensemble average and we may conclude 


._ |Bl’o*Aw 
ex on <G,,(x, x ‘): GS .a(x, x )>. 


We now drop the uninteresting pre-factors and take a spatial Fourier transform to 
define a quantity which shall play the role of a diffuse-field energy density. 


E(Q, A) = | d3 xe”) G(x, x’): G*,o(x,x’)). (4.5) 


E should, by reason of statistical homogeneity, be independent of source position x’. 

On sufficiently long time scales (Q « ac) and long length scales |A| « a, E according 
to our main hypothesis, should obey a diffusion equation. It is therefore expected, for 
small Q and A, that E should have an Q and A dependence of the form 


constant 
= ——___ 4.6 

iN+ DA? (4.0) 
characteristic of the Fourier-domain solution of a heat equation (1.1), and where D 
is the desired diffusivity. 

The constant in the numerator of (4.6) may be evaluated by means of a Ward 
identity (ZIMAN, 1979 ; ECONOMOU, 1983). Note that 


E(Q,0) = (| atacucw, x): G*, o(x, ‘) 


x’ = x” 


If G is expanded in terms of its unknown, but real, normal modes u” each with 
natural frequency o,,, 
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G(X’, X) = dui(x MID To. . ip? aie 


with orthonormality condition 
l d? xu"(x) *U"(X) = dy 


then 


ud} (x’ ud (x)ad}' (x)ud" (X") ) 
E82, 0) = (ry . *T(@+ ie)? —@2][(@—ie+Q)? — 02) -¥ 


iN ui (x’)u x") AE a a, ) 
= (25. (w+ie)’—a@, (@—ie+Q)*—a, |/ - 


5 (Goi, x”) —G* +Oii (x’, x") >» =x" (4.9) 


The A = 0 value of the mean square response E has been found to be related to a 
simple mean response. The numerator in the last line of (4.9) can be recognized, as 
Q = 0, as the imaginary part of the trace of <G), divided by the total system volume. 
As such it is related to the modal density in the effective medium (ZIMAN, 1979). As 
the modal density will be perturbed from its bare value by terms only of order v’, 
that modal density can be well approximated here by its bare value. More formally, 
(4.9) may be evaluated by replacing <G(x)> with the inverse Fourier transform (3.13) 
of the expression derived for <G(p)> in the preceding section. 


l l 


E(Q, (0) = 202+Q? (2x)? d’p{g: —GJo+0 + 2g), — 2940 





i ‘i 202+; -—<¢, 
~~ 2@Q 2n? p* dp (w? —p?c7 +0,)(w* —p’c7 +07) 


(4.10) 





20Q+607-—6r 
(w* —p?c7-+o7)(@* —p’ cz +47) 


+2 


where 7 has been neglected in comparison with 20Q. We may further neglect 2@Q 
relative to o*—o = —4iwac; that is, Q, the rate scale for study of diffusive field 
evolution, is slow compared to the rate of multiple scattering. We also approximate 
o <« w’, which is valid at all frequencies if the modulus fluctuation scale v is small. 


Equation (4.10) is then easily evaluated ; 
l 
The value E has been found to be inversely proportional to iQ, corresponding to a 


step function temporal behavior of the total energy. The relative contribution of 
longitudinal and transverse waves in (4.11) is consistent with the standard mode 
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counting methods (WEAVER, 1982) used in the conjecture (1.3). Equation (4.11) can 
serve as a check on any theory of mean square responses. 

Comparison of (4.11) and (4.6) allows one to conclude with an expression for the 
leading behavior of the diffusion propagator 


E(O,A) = 7 [I/c} +2/cH] (i+ DA”), (4.12) 


Of importance equal to that of the energy density is the energy flux. That flux may 
be written in terms similar to those used for the energy by considering minus the 
product of material velocity and stress: 


Ik sail 3 G;,(x, x’, 1)G, j, (x, x’, t) Cikyy (X). (4.13) 


The energy flux is summed over source polarization as was the energy. The modulus 
C(x) may be taken to be the Voigt average modulus C° to leading order in the 
magnitude of the fluctuations. 

Energy conservation takes the form 


é= —0, j, +source, (4.14) 


which may be corroborated by appeal to the governing equations (2.1). As above we 
construct a temporal Fourier transform of the energy flux associated with the band- 
passed response. 


x e110" BC) B¥(’)( — 1) Gui (X, X GS, (X, X’) (4.15) 


1@ 
= (2n) ikuv | B| 7A@<G.i)(X, xX’) Go +.avj,u(X x’)>. 


We drop the same prefactors as was done following (4.5) and take a spatial Fourier 
transform to define 


J, (Q, A) ess [ars "aie Gai (X, X’)GS +0, u(X, x’)»> (4.16) 


which shall serve as our definition of diffuse field energy flux. Energy conservation 
(4.14) then takes the form 


iIQE = iA, J, +source. (4.17) 


While the diffusivity D may be defined by means of (4.12), it will prove more 
convenient to define it by means of a constitutive-like relation between energy flux 
and energy gradient: 


J, = iDA,E (4.18) 


hypothesized to be valid at small Q and A, and hypothesized to follow from a study 
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of the multiple scattering. It is readily seen that (4.18) and (4.17) together imply (4.6) 
and therefore (4.12). Thus the two definitions of D are equivalent. 


5. MULTIPLE SCATTERING THEORY FOR MEAN SQUARE RESPONSES 


Both the energy and the flux may be formed by contractions on Green’s dyadic 
covariance 


(Gip(X, x )GRy, y’)> 


which is a dyadic generalization of a scalar quantity for which there exists a multiple- 
scattering formalism (FRISCH, 1968; ECONOMOU, 1983; KIRKPATRICK, 1985). The 
covariance is a fourth rank tensor function of four independent spatial variables. We 
again suppress the wm and w+Q subscripts, henceforth taking the asterisk to indicate 
a complex conjugate and also to indicate a subscript #+Q rather than w. It proves 
convenient to study the covariance in a spatial transform domain, hence a fourth rank 
tensor function of four wave vectors is defined. 


Hid? (p+q —q—p’) 


= (2n)- ‘T[ [fe xd3x’ dy d?y’ (Gig (x, XGHy, yen ta tay (5.1) 


The energy density may be expressed in terms of H by 


E= | [oraa'rs.atttsa608y (5.2) 


l 
(2n)° 
which follows from (5.1) and the definition of E. Similarly the energy flux may be 
expressed 


J, = 5 my? ctu | a qd ‘pipe t+ And sa? H48, 15 jp. (5.3) 


The Dyson equation governed the mean Green’s dyadic. The Bethe—Salpeter equa- 
tion (FRISCH, 1968) governs the covariance. In the Fourier domain it is 


a: AN = ‘T? a0” (p— a+ fa Ss it eta Pak Kissa ae ae (5.4) 


It is written in terms of the double-mean-field Green’s function 


Dip+a = (G,,(P)><GE(p+A)>, (5.5) 


where <G) is as given in the previous section. 

The Bethe—Salpeter equation is formally exact, but approximations are required 
for the intensity operator K. In what is often called the ladder approximation because 
of the shapes of the associated diagrams K is expanded to leading order in the 





Diffusivity of ultrasound in polycrystals 73 


perturbations, in this case, to leading order in v. The approximation is similar to the 
FOSA and a dyadic version of FRISCH’s (1968) equation (4.36). 


PakKeva © = H(p—Ss) pps, (pi + A))(s; +4 je. (5.6) 


The specifications (5.5) and (5.6) complete the prescription for the calculation of the 
covariance H. The Bethe—Salpeter equation (5.4) is however very complex and further 
simplifications and approximations are required before it becomes tractable. 

The first set of simplifications may be obtained by focussing on a different dependent 
variable, S, instead of H, where S is a contracted version of the source of the covariance 


and is defined by 


Sin(p, A) = fe qi,(T 8 ab as Ties a5p;, (5.7) 


mL 2-4 = - (G(P)> a2 <G*(P+A)> ni’ - (5.8) 


By performing suitable contractions upon (5.4), the Bethe—Salpeter equation becomes 
Sin(p, A) = Sunt [a's mM AmK e+ ajl s+a Sy (S, A). (5.9) 


S contains all information required for the present purposes. In particular, we find 
that the energy density and flux are given in terms of S by 


l 
E= Qn): [a piucT ie. Sin (p, A) (5.10) 


= Coakit a+ S?.(p, A). (5.11) 


In order to obtain the diffusivity D by means of (4.18) it will be sufficient to obtain 
E and J as expansions in powers of A through the linear terms. This will allow the 
derivation of an expression for D and also allow verification of the constraint (4.11) 
and the constraint that J vanish at zeroth order in A. To obtain this expansion it will 
be sufficient to work with expansions of T and K and S to similar, linear order. This 
will provide a further simplification of the Bethe—Salpeter equation. While doing this 
expansion it is wise to recognize that the system provides two different wavenumber 
scales, w/c and the attenuation «. In the limit « « w/c, a constraint invoked in previous 
sections, and valid throughout the entire frequency range if v is small, we recognize 
that our expansion through terms linear in A should retain terms like A/« ~ Awc/Imo 
in comparison to unity, but may neglect terms like Ac/@. 

The assumptions a « w/c and A « w/c allow the Bethe—Salpeter equation to be re- 
cast as an equation of radiative transfer. The further approximation A « « provides 
the diffusion limit. This two step procedure is the business of the next section. 
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6. RADIATIVE TRANSFER 


In the above three expressions the integrals over the magnitudes of the wave vectors 
are all highly singular, as can be seen from the form of I, Eqs (5.5) and (3.22) while 
recalling that the self energies o, and in particular, their imaginary parts, are small. 
This observation leads one to attempt to approximate the integral by the contribution 
from the vicinity of the poles. Such an approximation is equivalent to replacing the 
factor p*T with 


pil im+A sa R'(p, A) (8.5 — PP p) (im — PiPm)O(\p| | \co|/cr) 
+ R°(p, A)(P.PpPiPm)O(\p|—|@|/e.). (6.1) 


This approximation is valid only if Imo « w’, or equivalently, « « w/c. Note that the 
unit vectors in the direction of p+A have been approximated by p, in accordance 
with the assumed limit A « w/c. 

If we were to further seek the diffusion limit A «a and Q « a/wce we would 
approximate the energy propagators R by 
con | iwQ 


R’ (p, A) = aad imo? + iwA*pe,/Imo;+ o° | ; 


2c3Im mo, 


L ox |, | 
R*(p, A) = - . “— |i Imot + b@A: pe, /Imo,+.. | (6.2) 


We now recognize that, as the wavevectors p and s will have magnitudes of order 
«/c, the terms in A in (5.6) and (5.11) may be neglected. Use of (6.1) and the neglect 
of the As leads to great simplification of the Bethe—Salpeter equation. It is now defined 
on the unit sphere and all A dependence in the kernel is by means of the factors R. In 
direct notation it is 


S(p) = 1+ [oragntee, ‘gS (8e0/c,) RY (8) + pKsoc7-1-&:S(Se/c7)R’(8)]. (6.3) 


The A and Q dependencies in S and R have been suppressed. After the approximation 
(6.1), the expression for the energy density becomes 


E = (2n)~? | a°st:t'sGoye,)R" (8) +2 8:S(a/c,)R’ (8). (6.4) 
The expression for the energy flux also reduces to an integral over the unit sphere. 
= (2n)~ ys Catt {a *S[(S [Cy )S5gSiSm Sh (Se/c,) R* (S) 


+ (&/¢1) (Sap — $,5) (Sim — S:8) Sm (80/7) R"(8)]. (6.5) 


At this point the equations can be made to admit of a simple physical interpretation : 
S(p) represents a source of radiant intensity in the p direction. It includes both primary 
and secondary contributions; R represents a propagation of that intensity. Equation 
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(6.3) indicates that the radiant intensity has two sources: I represents the primary 
source, the original impulsive force. The integral sums the secondary contributions 
due to scattering of intensity into the p direction from all other directions §. Thus 
(6.3) is an equation of radiative transfer. 

There is a considerable body of literature on the subject of electromagnetic radiative 
transfer (CHANDRASEKHAR, 1960; SOBOLOV, 1963 ; ISHIMARU, 1978). The bulk of it is 
confined to a parameter range of little relevance for the present purposes. In that 
literature the effects of radiation polarization (the tensor character of the equation) 
are usually neglected, as is the time dependence Q. Any polarization which that 
literature does consider is, furthermore, confined to transversely polarizable waves 
and does not include the longitudinal modes implicit in (6.3—6.5). On the other hand, 
that literature is generally not confined to the diffusion limit A « «, nor to infinite 
statistically homogeneous media. It considers secondary sources and sinks other than 
those due to simple scattering and allows frequencies, w, to change upon scattering. 
Thus our present problem is both somewhat more complex and substantially less 
complex than the usual problems considered in the literature of radiative transfer. 

Equation (6.3) is a set of nine coupled integral equations for the nine components 
of the tensor S. As such it appears to remain somewhat intractable. However, our 
purposes are such that only two of those components will prove to be relevant. One 


might write, in general 


3 


S(p,A)= >) p*p'Sh, (6.6) 


AM = 1 


where p' = p and p’ and p’ form an orthonormal triad. Substitution of this form for 
S into (6.3)—(6.5) shows that the four components S$} Sj, Sj and S} do not contribute 
to J or E or to the secondary source of radiation intensity. They may or may not 
vanish, but they are in any case decoupled and irrelevant. The irrelevancy of these 
four components may be summarized by requiring that S(p) have as an eigenvector 
p itself on both the right and the left. The decoupling is due to the vanishing of the 
cross terms g’g*’ and g’g** in the approximation (6.1). That vanishing is in turn 
traceable to the different speeds with which transverse and longitudinal waves propa- 
gate. There is no coherent propagation of interference between waves of different 
velocity. Transverse modes can, however, interfere coherently and S} and S$} do not 
in general decouple from the radiative transfer equations. It is at first somewhat 
surprising that one would need four components to describe the intensity of a trans- 
verse wave. One should, though, recognize that an elliptically polarized transverse 
wave requires specification of four independent quantities: the two principal inten- 
sities, the orientation of the ellipse and also the relative phase. These four quantities are 
equivalent to the four “‘Stokes”’ parameters used in the treatment of electromagnetic 
polarization in the theory of radiative transfer (CHANDRASEKHAR, 1960). Thus the 
present system may be said to require a fifth Stokes parameter as well, to represent 
the intensity of the longitudinal waves. 

S may be reduced further, from five independent components to three if it is recog- 
nized that in an infinite statistically homogeneous and isotropic medium any second 
rank tensor function of two vectors p and A must be a linear combination of the five 
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dyadics pp, I, AA, pA and Ap. The condition that S also take p as a right and a left 
eigenvector then reduces this to three dyadics. The most convenient representation is 


of the form 


A(p, A) ... B(p, A 
RE PP) R! 


where A, B and C are isotropic scalar functions of p and A, and therefore of p*, p: A 
and A? (as well as Q). Since, as discussed above, the product RS represents radiation 
intensity we see that A, B and C represent radiation intensity; B and C represent 
transverse wave intensities ; A represents longitudinal wave intensity. 

The energy density and energy flux take very simple forms in terms of A, B and C: 


sa Ci 
Sip, A) = ip + (pp) Alta —pp) - a) —P 


E [a°8(4G0(e,)+2B160)e7)+ (8° A)}CBolen) (6.8) 


~ (2n)3 


J | d?§ §{c, A(8@/c,) + 2c7B(Sa@/cr) + c7(A* —§* A*)C(Sm/cr)}, (6.9) 


~ (2n)? 


where the A and Q dependences in the Stokes-like parameters A, B and C have been 


suppressed. 
When the form (6.7) is employed on the right-hand side of the equation, (6.3), of 


radiative transfer, it becomes, in direct notation, 


S(p) = I+ | d? $9 (p—Se/c, a? /cZ E28 A(S/c,) 


+ [°HP—So/ep.0*/eFE "BA —§)BSa/er) + f-$):4CSa/cr)]. (6.10) 


Equation (6.10) is equivalent to three coupled scalar integral equations defined on the 
unit sphere. They are obtained by contracting (6.10) on the three different dyadics pp, 
I—pp and (A-(I—pp))(A-(I—pp)) respectively. The first two are, with terms of 
explicit order A’ neglected, 


A(poo/c,)/R’ = 1+ | d?5n' (0,,)e0* cf (RRS...) A(S/c,) 
w* ans ‘ 
+ | d?5n'"(0,,) = BRU 3): B(S@/cr)] (6.11) 
CLCT 
and 


a » a* inet cate a ek 
2B(pa/cr)/R’ = 2+ |a°sn"(,) al ((1—§) Bs-:::A (Seo/c,)] 
LCT 


+ [asm **@..) = (1 —8)8 (1 —$)::B(Sw/cr)]. (6.12) 
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The third equation is not required for the present purposes and is therefore not 
included here. The angle @,, is the angle between the p and § directions and the n(0) 
were defined in (3.33). Except for the factors of A and B in the above integrands, 
these integrals are equivalent to those of equations (3.31)—(3.35). 

In the diffusion limit, where A « «, we seek a regular expansion of the Stokes 
parameters A and B in powers of A/a. As discussed above, retention of A dependence 
through the linear terms will suffice. Hence we wriie 


A(po/c,, A) = Ay+A,(A*p) (6.13) 
and 
B(po/c,, A) = By +B,(A°p), (6.14) 


where the coefficients Ao, A;, By and B, are independent of A. As they are isotropic 
scalars they must also be independent of p. They may, however, retain Q dependence. 
The energy propagators R must also be expanded ; 


R’(p) = po +pi(A:p) 
R*(p) = po +pi(A-p) (6.15) 


whence, upon reference to (6.2) we identify 


a (1- =a () iwQ 
Po 2¢3Imat Imo% Acta, 2¢,a,/)’ 


= (1- )- ed (1- iaQ 6.16 
po 2¢3Ima™ Imo’)  4cta, 2cra,)’ ney 











and 


iw?n in 


- o- ia 
2 —_ ’ 
2cz(Ima,)? — 8c}a? 





‘ iw?n in 
Daas = . 
2c7#(Imo,y)? — 8c4a2. 





p (6.17) 


The diffuse field energy density (6.8) is given through first order in A in terms of 
the expansion coefficients of the field by 


I 
E = 573 (Ao + 2B). (6.18) 


The energy flux (6.9) is given through order A’ by 
J = (A/6n’)[A ,c, +2B icy]. (6.19) 


By reference to (4.18) the diffusivity is then given by 
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i Aye, +2Byer 
3 Agt2B, — 


D=- (6.20) 


Were we to seek the diffusivity by comparison of E with the expected form of the 
diffusion propagator (4.12) instead of by means of comparison of the flux with (4.18), 
E would have been required through second order in A. Correspondingly, A and B 
would have been required through that same order, and C to zeroth order. By 
considering energy flux and (4.18) instead of (4.12) an energy conservation constraint 
(4.17) is implicitly being employed to aid in the solution to the equations of radiative 
transfer. 

The integrals in (6.11) and (6.12) over solid angle § will be carried out using 
conventional spherical coordinates 6,, and #,, referred to a polar axis in the direction 
p, and to an arbitrary zero of longitude. The factor §*A = Acos@,, appears in the 
integrand and must be re-expressed in terms of the chosen spherical coordinates. This 
expression is available from spherical trigonometry, or alternatively and in greater 
generality, from the addition theorem of spherical harmonics. It is 


s-A = Acos 6@,, 
= A[cos 0, cos 8,4 +sin 0,, sin 8,4 cos [b,, — Pap] 
= (p- A) cos @,,+azimuthal terms. (6.21) 


The azimuthal terms will vanish upon integration over @,, in any integrand such as 
those of (6.11) and (6.12) which are otherwise independent of ¢,,. This was to be 
expected. There are no terms on the left-hand side depending on sin @,,, hence such 
terms should not appear on the right-hand side. 

The expressions (6.13) through (6.15) are now inserted in the integral equations 
(6.11) and (6.12). Like powers of p-A are collected and the integrations carried out 
in terms of the definitions (3.35 ff) of the attenuations and their angle-weighted primed 
versions. From the A-independent terms one deduces 


Ay = pot (4pbc7/2) [4110140 + 2a7,¢7Bo] 


2Bo = 2p$ + (4phc7/n)[a, 70, A 9 + 2a77¢7Bo]. 
From the terms linear in p- A one finds 
A, = pt +(4picr/n) [41161 A0 + 2471¢7Bo] + (4poe7/m) [acer A 1 +20%7,¢7B)] 
and 


2B, = 2p) + (4picr/n)[ar7c, A 0 + 2477¢7Bo) 
v (4picr/n) [a 7c, A | +2a77crB,]. (6.23) 


To leading order in Q the solution to (6.22) for the lower order terms is 
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n | n | 
2iQ c}. 


Ay By (6.24) 


2iQ ¢}’ 


This result together with (6.18) corroborates the expected total energy constraint 
(4.11). It serves as a valuable verification of the extensive algebra of the present section 
and also of the consistency of the ladder approximation. 

At the next order in A and again at leading order in Q, we find 


, , 3 
™ (Ap — App + Ap7)/C7 


A => , , , , 
1 40 (4, — OH), )(47 — Ayr) — Ap 7X7, 


, 1 (%, — Op, +07,)/e7 
‘ia (a, — Op, )(Gp—A py) — pay 
Hence, using (6.20), 


(a7 — Apr + a1 7)/Ci +2(a, — a), +4 7,)/c7 


6 ( l 2 ) (4, —O1,)(A7—Ayy) — Ap 7X7, 


a’ 

This is the desired result. The diffusivity of the multiple-scattered wave field is 
expressed in terms of the wave speeds, the attenuations and the scattering-angle- 
weighted attenuations. It is valid for statistically isotropic and homogeneous media 
in the limit of weak attenuation per inverse wavenumber, « « w/c, and equivalently, 
weak fluctuations, v, in moduli. It has been established for frequencies below the very 
high frequency, geometric optics regime, that is, within the “Born” approximation. 

In the Rayleigh limit where all primed attenuations vanish, the expression for the 
diffusivity simplifies : 


D o=(2 re. sr \i(5 +3) (6.27) 
— c} 6a, cp 6ar a Fi J 


which is precisely that diffusivity argued as plausible in the first section. It is striking 
that the naive rule of mixtures used there was correct. 


7. NUMERICAL RESULTS AND DISCUSSION 


Explicit expressions for the attenuations, and therefore for the diffusivity, can be 
obtained only upon specification of a two point correlation function n(r). The choice 
(2.5) has the virtue of simplicity while adhering to the requirement that 9 vanish at 
large r and reach v’ at r = 0. That specification has a single length scale B~'. In order 
to make comparisons with real materials one must either determine the actual n(r) of 
the specimen at hand (STANKE, 1986) (and by implication verify the decoupling 
assumed in (2.4)) or make a best estimate for an effective 8. At low frequencies it is 
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well understood that the key microstructural quantity is mean grain volume. Hence 
one identifies an effective B by the equation 


V = favre = (2n)°4(0)/v? = 82/B°, (7.1) 


where V is the mean grain volume. At frequencies of order w ~ fc or higher, though, 
it is probably not the grain volume which is critical, but rather the short distance 
behavior of 7, or the volume density of grain boundary. This leads to the identification 
of # with half that density. Such a Pf is something like half the B determined from 
(7.1), depending on the assumed grain shapes. This uncertainty in the choice for an 
effective 6 underlines the naivete of the simple form (2.5). 

We shall, though, in order to generate explicit expressions for the diffusivity, make 
the choice (2.5). It should be kept in mind, however, that in any comparison between 
the results to be presented and diffusivities in real materials, the best choice of an 
effective f is uncertain, and should, probably, even vary somewhat with frequency. 
With choice (2.5) the spatial Fourier transform of is obtained from (3.10). 


i(q) = (v°B/n*)(B +4) *. (7.2) 


The n(@) are then found from Eqs (3.33). 
n“(0) = (vin?) 


n'™ (0) = (vin?) 








w 


n'*(0) = n°"(0) = (v°B/n’) | B+ oi aa 


CT CL 
which in turn allows the attenuations, (3.35) to be written as 


B y? (*1 
ue = 1950 of Xt | 


> B y2 ; *1 
Orr = 2100 ca? J 


| du(3+p’)/[1+2x7(1—w)? 





du(24+ 3y? + w*)/[1+2x7(1—p))’, 
1 


er 
ur = F950 ¢ kext | AMIS + 6u?— mV + x7 +x7—2xrxrW))’, 


tr, = a, 7c7/2c; (7.4) 


where x, = w/c, B and x; = w/c;f > x, are dimensionless measures of frequency and 
u has been used for cos 8. The primed attenuations are given by expressions identical 
to those above, but with one extra factor of y in each integrand. 
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The integrals are easily done in the Rayleigh limit where x, < x; « 1. The primed 
attenuations vanish in comparison to the unprimed attenuations; the unprimed 
attenuations are 


app" = (8/375) By? x? /cz, 
app = (9/375) py? xz/cr, 
app 8 = (12/375)pv?xixz/crcr, 
opeyeeh = (6/375) Bv?x7x7/crcz (7.5) 
and the attenuations proper are (3.34) 
ORayteigh = (1/375) B(v*/c7)x7(8 + 12K°), 
ORayteigh = (1/375) B(v*/cr)x7(9+6K ~*), (7.6) 


where K is the ratio of Voigt wave speeds; K =c,/c;. These attenuations are in 
agreement with expression found in the literature (STANKE and KINo, 1984, Eqs 24 
and 25). 

The Rayleigh limit diffusivity is then given by (7.6) and (6.27) 


as OE a Neg htalie Vii (; “) 
Drayleigh = 6B" (v /cr) (2+ K~*)(1+3K~°) 3 + 4 ° (7.7) 


The result is.nearly insensitive to K in K’s physical range. For K about 2, Daayieigh 1S 
given fairly accurately by 


Drayleigh ~ 11.2(¢r/B)xr *(v*/er)” : (7.8) 


It is noteworthy that the transverse and longitudinal contributions to the energy flux 
are roughly equal. 

Outside the Rayleigh limit the integrals (7.4) are more difficult. They are all of the 
form 


a 1 yu" 
tw | on aig ni 


forn =0,1...5 and |b| < |a|. Defining ¢ = a/b the integrals are given by 














3) o+l 

I, = pb igre ett 
_2) 26’ p+ 

I,=b gin ~2#lont +t 42h, 
2) 2¢° 21.0?%+! 

I,=b ae 6 oat *t 40h. 
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In terms of ¢, = 1+1/2x; the longitudinal attenuations are 


pw pas o. +1 


o 6? + 38/3 
4200 ci] 2-1 $,-1*°ret 


Are —4,(3+¢7) log 


and 


po. +l 


—(9+ 18; + 51) log | 
ee 


/ 
Ars 


— B v?| 26,[3+ 7)’ 
~~ 4200 c3 


86} +76¢,/3 |. 
ob; —1 7 pi + py, 


(7.11) 


The transverse attenuations are, in terms of @; = 1+ 1/2x;, 


B vw | 2B eo br+l 


= : 663+ 20/3 
“TT * 8400 c3 p7—1 =i Ort 


—27(3 + 23) log 


and 


prt 


»2 4 
B 4 | 2PrA SIO) (24449934 561) 08 
ae 


= +8434 40¢,/3 |. 
“TT ~ 3400 cé or—1 sill | 


(7.12) 


The longitudinal-to-transverse mode conversion attenuations are, in terms of 
, 9 
Pr. = (1+ x7+7,)/2X,X7, 


ae es | 25 +h — 91) 
“T ~ 4200 cc} o2, —1 


ori v I 


apy — 667, +34/3 
) 


—4o7,(3- pin) log 


and 


B v 74x (1 5+ 6b 71 ni pir) 


o7,—1 


— (15+ 187, —S$7,) log 


“LT 4200 c.c3 


Prti 3 Pr 
$7, —1 oon +85 | (7.13) 


The other mode conversion attenuations may be found from 


tr, = a&,7(c7/2c7) 
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aggregate of cubic crystallites of weak anisotropy. The low frequency Rayleigh regime is characterized by 

an inverse scaling with the fourth power of frequency. The high frequency stochastic regime is characterized 

by an inverse logarithmic dependence. The four data points from the experimental work of Guo et ai. 
(1985) are also shown. 


Mr, = Hr(c7/2c7). (7.14) 


As x—oo the resulting expressions for the attenuations «” = «,,+a,; and 
a’ = a&;rp+4,7, approach the “stochastic asymptote” and agree with expressions given 
by STANKE and KINo (1984) in their Eqs (33) and (34). 

The diffusivity is obtained upon substitution of (7.11) through (7.14) into (6.26). 
The result is plotted in Fig. | in scaled dimensionless units of diffusivity against 
dimensionless frequency for a choice of Poisson ratio 0.287 corresponding to a wave 
speed ratio K = 1.827 appropriate to untextured polycrystalline iron. The diffusivity 
is given in units of 4200c?/v’B; frequency is measured by x; = w/c7f. The inverse 
fourth power dependence is apparent in the Rayleigh regime x; « 1. An inverse 
logarithmic dependence is seen in the opposite limit. The diffusivity must certainly 
become frequency independent in the very large x, geometrical optics regime. Thus 
the observed inverse logarithmic (and therefore weak) dependence closely, but not 
precisely, captures the actual very large x limit. 

The approximate point at which geometrical optics becomes significant and the 
“Born” approximation fails can be determined from the attenuation curves published 
by Stanke and Kino. For the case of iron, with a moderate value of v that point is 
found to be at x, = 4.5; for aluminum, with a small value of v one find that geo- 
metrical optics obtains for x; > 20. One speculates that the diffusivities should become 
frequency independent for frequencies in excess of these limiting points. It is note- 
worthy that the calculated diffusivities are nearly frequency independent there anyway, 
in spite of being based on attenuations calculated using the Born approximation. 

The analogous situation for the attenuations was very different and should be 
contrasted with that for the diffusivity. The stochastic asymptote has attenuations 
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increasing with the square of frequency while the geometrical optics very high fre- 
quency limit has frequency independent attenuation. Thus the stochastic asymptote 
is a very bad approximation in the geometrical optics limit if the interest is in 
attenuation. The error is apparently much less severe in the case of diffusivity. 

The four diffusivities measured by Guo et al. (1985) are also plotted, as isolated 
points. In order to do so, their measurements had to be nondimensionalized. A shear 
wave speed of c; = 3.23 mm us _' was chosen corresponding to untextured iron and 
v/c; was set to — 1.66 in accordance with the moduli of single crystal iron. A value 
of B was chosen equal to 10 mm~'. The one-parameter fit represented by our free 
choice of a value of f is seen to agree fairly well with the limited laboratory data. 

The value B~' x 100 microns appears, from the optical micrograph published by 
Guo et al., to correspond to the radii of the regions of similar shading. It is reasonable 
to conjecture that the grains clump in groups with similar crystallographic orientations 
and thus with similar shadings in the micrograph. The ultrasonics would of course be 
most sensitive to the large angle grain boundaries between the clumps, and hence 
respond to an effective B~' of the order of 100 microns. It is, unfortunately, not 
possible to make detailed comparisons without more information on the actual sample 
microstructure and alsc on error bounds on the diffusivities. The agreement is, 
however, promising. 


8. CONCLUSIONS 


The present theory has presented a formula for the diffusivity of a multiple-scattered 
diffuse wavefield in an untextured cubic polycrystal. It appears to agree with lab- 
oratory measurements. The theory is valid in the limit of weakly anisotropic crystallites 
and is derived within the context of the “Born” approximation which thus confines 
confident application to the so-called Rayleigh and “‘stochastic’’ regimes where the 
frequency lies below the geometrical optics limit. The diffusivity was derived in a two- 
step process. Radiative transfer equations were obtained from the multiple scattering 
equations in the ladder approximation by invoking a weak attenuation and long- 
length scale assumption. The diffusion limit was then obtained from the radiative 
transfer equations by considering still longer length scales. 

There is room for further study on a variety of issues. Perhaps the most obvious 
need is for definitive laboratory studies comparing the diffusivities in known micro- 
structures with the predictions Of a theory such as the present one. Inasmuch as it is 
as yet unclear how small the numerical quantity |v/c7| must be for the theory to be 
valid, it would be desirable to study diffusivities in a range of materials with different 
values for that parameter. Such experiments would also tend towards the development 
of efficient techniques for the measurements of diffusivities, and by implication, local 
internal friction as well. The possible consequences for NDE are exciting. 

It is conceivable that the formula for diffusivity could have been derived from a 
starting point within the theory of radiative transfer without any direct reference to 
the formalism of multiple scattering, with a consequent saving of effort, and with a 
consequent loss of grounding in the fundamental mechanics. Diffusivity in systems 
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with non-trivial geometries or with statistical anisotropy might be best treated in that 
manner. 

The present calculation has been carried out within the “Born” approximation and 
is therefore not safely applied to the geometrical optics very high frequency regime. 
While attenuations have been calculated in that regime (STANKE and KINo, 1984; 
KARAL and KELLER, 1964) they are of little practical relevance, being of order f itself 
and too large for meaningful laboratory measurement. Laboratory diffusivity in such 
a regime, though, is not inaccessible; indeed the geometric regime has even been 
termed the “diffusion regime”’ by some authors. Such an appellation is somewhat of 
a misnomer in that diffusion is operational at all frequencies on sufficiently long length 
scales. It is appropriate though in that no other ultrasonic microstructural assessment 
method is likely to be feasible there. Thus it is unfortunate that the present calculation 
cannot be applied there. There is a need therefore for a diffusivity calculation in the 
geometrical optics limit. 

In the limit of strong crystallite anisotropy, FOSA and the ladder approximation 
will fail, as will the radiative transfer approximation « « w/c. The condition « ~ w/c 
is precisely the Ioffe—Regel criterion for the absence of diffusion and the onset of 
Anderson localization (ZIMAN, 1979; ECONOMOU, 1983; KIRKPATRICK, 1985; 
WEAVER, 1989a). While an analysis of the high crystallite anisotropy regime is pre- 
sumably beyond any exact treatment, experiments there would not be difficult and 
should certainly be attempted. 

With a theory now available it is possible to make the case that ultrasonic diffusivity 
measurements should take their place with attenuation and backscatter measurements 
for purposes of nondestructive characterization. While the technical aspects of making 
diffusivity measurements are non-trivial, and will be discussed in another com- 
munication, it seems that the severe systematic errors which can sometimes plague 
attenuation measurements do not trouble diffusivity measurements. The ease with 
which quantitative assessment of microscale length was effected at the end of the 
preceding section is encouraging. The diffusivity measurements also promise to deliver 
assessments of the local internal friction, something which it is difficult to imagine 
measuring in any other non-destructive way and which could conceivably be an 
interesting nondestructive materials evaluation parameter in its own right. 
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ABSTRACT 


EXPERIMENTAL Strain distributions are determined very near the crack tip in Fe-3wt.%Si single crystals. 
Both in situ stereoimaging and electron channeling techniques give reasonably reproducible distributions. 
By growing fatigue cracks on a {100} cleavage plane, the singularity strengths have been determined for 
both growing and stationary cracks under relatively plane stress and plane strain conditions. This has 
allowed a comparison to existing theoretical models. It is shown that the HRR singularity (Hutchinson, 
Rice and Rosengren, 1968) for stationary cracks is very good to within | wm of the crack tip and a 
hardening mode! for the growing crack (GAO and HWANG, Advances in Fracture Research, edited by D. 
FRANCOIS, 5th Int. Conf. on Fracture, Cannes, France, 2, 669, 1981) is surprisingly good. Other issues 
such as fracture criteria are discussed since strains greater than unity were measured at the crack tip in this 
relatively brittle material. 


INTRODUCTION 


THE IMPORTANCE of knowing the distribution of strain within the near tip region of a 
crack is substantial and considerable effort has been expended to determine what the 
strain functions are. Most of the work on determining strain distributions has been 
analytical, with only a few experimental efforts made to substantiate these analyses. 
It is not known how well continuum solutions are likely to apply in a region which 
may be dominated by dislocation emission from the crack tip (RICE and THOMSON, 
1974). For example, a crack tip with a plastic zone of 1000 4m would have a computed 
crack tip opening displacement of 2.5 um, which would require a large number of 
dislocations to be emitted from the crack tip. These dimensions generally agree with 
the simulation work of ZHAO et al. (1985) using dislocation emission and equilibrium 
criteria for dislocations distributed on two symmetric slip planes about the crack tip. 
For these crack tip conditions, the simulations indicated that dislocations would be 
spaced on the order of 0.1 wm apart near the crack tip, and that for distances greater 
than | ym, the continuum solutions should apply. This large number of shielding 
dislocations changes the stresses near the crack tip, and thus the strain distribution. 
Just how close to the crack tip this redistribution of strain occurs is likely to be 
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important from the standpoint of both ductile and brittle fracture mechanisms which 
depend on a critical distance parameter (RITCHIE and THOMPSON, 1985). Molecular 
dynamics solutions to near crack tip events (DAw and Baskes, 1983) also lie within 
the region which may not be applicable to continuum mechanics solutions. Defining 
just how large a region might be appropriate to continuum solutions would help to 
limit the number of atoms in discrete atom simulations of crack tip events. 

Most of the analytical and numerical solutions which exist for growing cracks are 
for elastic, ideally plastic solids. For strain hardening solids two studies have been 
made. A mode I asymptotic analysis in plane strain has been made (AMAZIGO and 
HUTCHINSON, 1977) for a material with a bilinear stress—strain curve, and a similar 
analysis has been made for a solid having the flow stress simulated by a Ramberg 
Osgood equation (GAO and HWANG, 1981). However, the relevance of the first of 
these is reduced for the growing crack (LAM and MCMEEKING, 1984) because of 
neglect of plastic flow along the crack flank, and the second has been difficult to 
corroborate (HUTCHINSON, 1987). Even if the continuum solutions for stationary or 
growing cracks are appropriate to the near crack tip region, there is little experimental 
confirmation that they are numerically correct. 

Early experimental work to determine the strain distribution near crack tips was 
hampered by insufficient spatial resolution to obtain strains within about 100 um of 
the crack tip. Photoelastic coatings (GERBERICH and SwWEDLOw, 1964) and Moire 
methods (Liu and IINo, 1969; UNDERWoobD et a/., 1971) are examples of these studies. 
Later, higher spatial resolution techniques were developed, using selected area electron 
channeling (DAVIDSON and LANKFORD, 1981) and stereoimaging (DAVIDSON, 1979). 
The latter technique was used to determine strains at the crack tip, while the former 
was limited by the severity of deformation at the tip. While the stereoimaging technique 
allows information to be obtained at the crack tip, it is limited to doing so only at the 
surface. Conversely, the selected area electron channeling technique can also be used 
to obtain information on the interior of the material (DAVIDSON and LANKFORD, 
1980), either by removing material from the surface of the specimen, or by making 
measurements below the fracture surface after passage of the crack. 

Only for low carbon steel (DAVIDSON and LANKFORD, 1981) have both the elec- 
tron channeling and stereoimaging techniques been used on the same material. A 
distribution of subgrains was found near the crack plane when the material was 
imaged with backscattered electrons on a section about 100 um in from the surface, 
or deeper. The sizes of these subgrains were measured along several directions from 
the crack tip, and it was possible to infer the stress distribution around the crack tip 
using the relation between fatigue formed subgrain size and stress (LAWRENCE and 
JONES, 1970). It was not possible to obtain strains accurately from channeling contrast, 
since the strain distribution is related to subgrain misorientation. 

The stereoimaging technique was used also to determine strains near fatigue crack 
tips in this material (DAvIDSON and LANKFORD, 1981). Strain distribution in the 
direction of crack growth was found to be fit best by the function 


( | 
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where Az is the fatigue strain, r = distance from the crack tip, and A, m and B are 
constants used to obtain a best fit to the data. Such a formulation also fits the strains 
at @ = +72/2 from the crack plane. This function fits the strain data marginally better 
than the relation 


Ae(r) = A/(r+ By)". (2) 


Two representative sets of data are shown in Fig. | for a 0.05 wt.%C low alloy 
steel tested in vacuum (DAVIDSON and LANKFORD, 1981). The dotted curve is similar 
to (1) and is only slightly displaced compared to the best statistical fit to the open 
circles. The second set of data (closed circles) was taken at a slightly reduced load 
ratio but the same AK. 

Strains ahead of crack tips in 7075-T651 and 7091 aluminum alloys, Ti-6Al-4V 
titanium alloy and 304 stainless steel have also been fit with (1) and (2), and it was 
found that (1) again fits the data better than (2) (DAvipson, 1986). In addition, the 
cyclic stress-strain equation for low-carbon steel was used to compare stresses derived 
by stereoimaging on the surface of the specimen with stresses derived from subgrains 
imaged 100 «xm below the surface, and good agreement was found (DAVIDSON, 1988). 

The experimental studies which have been performed to date are listed in Table 1, 
which also gives a representative list of the analytical and numerical studies. Clearly, 
more experimental information is required to examine the validity of the analytical 
and numerical work, and this is especially true for tensile-loaded cracks in plane 
strain. 

The only known experimental work which approaches plane strain for a stationary, 
tensile crack is the study of a Charpy V-notch specimen by (WILSHAW, 1966) and a 


recrystallization study of both stationary and growing cracks by HORNBOGEN and 
MINUTH (1981). Although the latter investigators measured strains in the center 
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Fic. 1. Strain distributions ahead of two fatigue cracks in a ferritic, low-carbon steel (¢, = 218 MPa) 

grown at AK = II MPa ,/m, R = 0.16 in laboratory air (DAVIDSON and LANKFoRD, 1981). Strains were 

determined by stereoimaging. The dashed line is a fit to Eq. (1) and the solid line is that of Eq. (7), modified 
for total strain. 





TABLE |. Analytical, numerical and experimental studies of crack-tip strain distributions 





Elastic, ideally 
plastic (e—p) 


Stationary — 
cracks With strain 
hardening (e—p”") 


Elastic, ideally 
plastic (e—p) 


Growing — 
cracks With strain 
hardening (e—p”") 


Analytical 


Plane stress 
Hult and McClintock 
(1957) 
Dugdale (1960) 
Hahn and Rosenfield 
(1965) 


Plane strain 
Green and Hundy 
(1956) 

Ewing and Hill 
(1967) 


Hutchinson (1968), Rice (1968), 
Rice and Rosengren (1968) 


McClintock (1958) 

Chitaley and 

McClintock (1971) 
Mode III 


(none) 


Slepyan (1974) 
Gao (1980) 

Rice, Drugan and 
Sham (1980) 


Amazigo and 
Hutchinson (1977) 
Gao and Hwang (1981) 
(no others) 


Numerical 


Rice and Johnson (1970) 
Levy et al. (1971) 

Rice and Tracey (1973) 
Lin and Chiang (1981) 


Kfouri and Miller (1976) 
Sorensen (1979) 


Dean and Hutchinson (1980) 
Lam and McMeeking (1984) 


Experimental 


Plane stress 


Gerberich and 
Swedlow (1964) 
Liu and lino 
(1969) 


Underwood ef al. 


(1971) 


Hornbogen and 
Minuth (1981) 


Plane strain 


Wilshaw (1966) 
Hornbogen and 
Minuth (1981) 
(no others) 


Davidson and 


Lankford (1980), 


(1981), (1983) 
Davidson (1987) 
(no others) 
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region of a 12.7 mm thick plate (plastic zone ~ 4 mm), insufficient information was 
presented to establish the strain distribution. Therefore, only very limited experimental 
crack tip strain distributions exist for stationary and growing cracks which can be 
used to examine the analytical and theoretical studies listed in Table 1. 

The focus of the work reported here was to examine crack-tip strain distributions 
at the surface and in the interior using the capabilities of both selected area electron 
channeling patterns (SACP) and stereoimaging. SACP are of limited spatial resol- 
ution, with the smallest beam size being about 5 um, and unuseable at large strains. 
However, SACP are derived from very near the surface—about 100 nm for Fe at 20 
keV; therefore, the spatial resolution in the depth direction is very high. During the 
determination of near-surface strain gradients, this high depth resolution may be used 
by depth profiling (sequentially electropolishing away the surface). Stereoimaging 
presents complementary capabilities to those of SACP in that this technique is capable 
of determining large strains with a high level of spatial resolution, typically 0.25 um, 
but only at the free surface. High purity single crystals of Fe-3wt.%Si were chosen 
for this study because of previous success achieved in growing planar fatigue cracks 
on a known plane, the (100), and in a known macroscopic direction, the [010]. This 
degree of control over crack growth was achieved by cycling the specimens in one 
atmosphere of hydrogen (KACZOROWSKI et a/., 1986). These experimental conditions 
resulted in the following benefits : (1) a minimization of crack growth in mode II and 
mode III, which could confuse the experimental results; (2) since the material was 
single crystal, inhomogeneous deformation would be caused solely by anisotropy 
rather than grain boundaries and second phase particles ; (3) cracks grown in hydrogen 
exhibit limited plastic deformation, which has two benefits, (a) SACP could be used 
nearer the crack tip, and (b) the effects of crack tip constraint are maximized. 

Anisotropic behavior is anticipated from single crystal, bcc Fe-3wt.%Si crystals. 
Slip occurs in the <111) directions on the {110}, {112} and {123} planes, for a total 
of 48 slip systems, with very little difference in critical resolved shear stress between 
them. Thus reasonably homogeneous plasticity with multiple slip is expected except 
at the earliest stages of yielding. 


2. EXPERIMENTAL MATERIALS AND TECHNIQUES 


In order to provide constrained plasticity, relatively large single crystal specimens 
4.8 mm thick were used. Single crystal material was relatively scarce, so specimens 
were prepared for crack growth on the (001) and in the [010], as depicted by the 
crystallography in Fig. 2. Specimens were wire cut from round single crystals, ground 
flat using 600 grit silicon carbide papers and then electropolished (KACZOROWSKI ef 
al., 1986). 

From earlier work (VEHOFF and NEUMANN, 1980) it was known that this crys- 
tallography of crack growth would produce asymmetric slip at the crack tip. To 
evaluate the severity of this asymmetry, an anisotropic elastic finite element analysis 
of the stresses on the 48 available slip systems was conducted. The crack tip region 
was divided into 16 sectors and evaluated above and below the crack plane. In the 
sectors from 45 to 135°, and from 225 to 315°, stresses on the five largest resolved 
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Fic. 2. Orientation of the test samples with a [001] load axis and a macroscopic [010] growth direction. 
Two of the 48 possible slip systems are shown. 


shear stresses were found to be about the same magnitude across the crack plane, and 
they varied only by about 10% within each sector at a fixed distance from the crack 
tip. In addition, each of these sectors had at least one representative system of each 
of the {110}, {112} and {123} types. Thus, some asymmetry was expected, but not 
much. 


2.1. Material properties 


Three material properties were required to accomplish a comparison of the exper- 
imental strain gradients to continuum models. The tensile modulus of elasticity in the 
direction of loading, Ej; 99, was 1.32 x 10° MPa. The yield strength and strain-hard- 
ening exponents were determined from compression testing solid cylinders with a 2:1 
height to diameter ratio. These 4.8 mm high samples were taken from the undeformed 
sections of fatigue specimens. Test parameters were ambient temperature with an 
applied plastic strain rate of é, = 0.005 s-'. This was chosen since it approximated 
the average rate of 0.006 s_' as deduced for crack-tip strain rates in the fatigue-crack 
propagation test (KACZOROWSKI ef al., 1986). Duplicate tests gave o,, = 296+ 
14 MPa and a constitutive flow equation given by 


Thow = 1000e)** (MPa). (3) 


While not accurate for plastic strains less than 0.04, Eq. (3) was good for 
0.04 < ¢, < 0.24. Note that the strain hardening exponent, 1 ~ 0.38, here, is 1/N 
where N is the exponent typically used for power-law hardening in the HUTCHINSON ef 
al. (1968) small-scale yielding formulations. Although these results are for monotonic 
compression testing, they are not largely different than results obtained by IKEDA 
(1981) for cyclic hardening of pure iron. Here ¢100> oriented Fe crystals gave a cyclic 
strain hardening exponent of approximately 0.33 at 150°C. It should be mentioned 
that pure iron does cyclically soften at such temperatures after a relatively few cycles. 
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This is associated with well-defined cell walls and no dislocations in the interior of the 
cells. Since no such cell walls were found (KACZOROWSKI e/ al., 1986), it is hypothesized 
that the cyclic softening stage was not reached in this Fe-3wt.%Si single crystal. 


2.2. Fracture mechanics test procedures 


Detailed experimental test procedures for growing cracks in hydrogen have been 
given elsewhere (KACZOROWSKI ef a/., 1986). Briefly, small compact round tension 
samples were loaded in tension—tension fatigue at cyclic frequencies ranging from 0.03 
to 3 Hz. The fracture surface produced, as shown in Fig. 3, was very flat. For a fatigue 
crack growing in hydrogen with Kj"** = 14.4 MPa ,/m, R = 0.1, the plastic zone size 
at the surface was estimated from surface deformation to be somewhat less than 
| mm. This is consistent with small scale yielding estimates of the plastic zone size, 
R,, given by (MCCLINTOCK and IRwIN, 1965) 


R, = Kj/no;, plane stress, (4), 
R, = Kj /3x0;, plane strain. (5) 


For plane stress, the estimate from (4), gave a calculated value of 0.75 mm for 
K, = 14.4 MPa ./m while from (5) it was 0.25 mm for plane strain. Thus, by this 
definition, experimental conditions approach plane stress at the free surface. After 
about 3 mm of crack growth at this AK (da/dN x 2x 10° m/cycle), two samples 
fractured and were subsequently examined using electron channeling. 

A matched sample, which was initially fatigue pre-cracked under these same con- 
ditions, was loaded 6 months later to K, = 13 MPa ,/m in a cyclic stage which fit 
within the scanning electron microscope. Stereoimaging was used to determine the 
strain distribution at this load level and then the specimen was loaded to K, = 25.2 
MPa ./m. After being loaded to this latter level, the specimen was found to have 
surface deformation at a distance exceeding 5 mm ahead of the crack tip. This large 
value indicated that a better theoretical estimate of plastic zone size would be obtained 
using the DUGDALE (1960) formulation 

= )- 1 (4), 


which gave a plastic zone size of 8.5 mm for K = 25.2 MPa ./m. Thus, the surface of 
the specimen under this tensile loading condition should be considered to be in a 
condition of plane stress. 


2.3. Electron channeling 


Selected area electron channeling patterns (SACP) were used to study the defor- 
mation associated with a growing fatigue crack. It was desirable to obtain information 
from both the surface and interior of the specimen. Two samples were examined, both 
from the fracture surfaces, as shown in Fig. 4. One sample was thoroughly studied, 
and the other only in the region near the surface. SACP were made from the locations 
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Fic. 4. Schematic of a fractured fatigue sample studied by electron channeling showing the locations of 
SACP relative to the starting notch and to the specimen surface. Surface was removed to a depth of 
0.8 mm on one side by cutting with a wire saw. 


shown in the figure. The sample studied the most had three regions defined by crack 
growth at AK = 14.4 MPa /m in | atmosphere of hydrogen with a loading rate of 
3, 0.3 and 0.03 Hz. 

The plasticity associated with fracture was examined for each of the locations 
between y; and y>, v2 and y;, and y; and y, in the figure. SACP were made at four 


different distances from the side of the specimen: x = 0.2, 0.5, | and 2 mm. The 
pattern taken from 0.2 mm was assumed to be within the plane stress region, while 
the other patterns were thought to be within the zone of plane strain. At each of these 
locations, SACP were made at three depths below the fracture surface. Material was 
removed from the fracture surface by electropolishing in a solution of 5% perchloric 
acid in acetic anhydride at 30°C. Use of 90 V and 30 mA in a jet polishing system 
which gives a circular depression approximately 2 mm in diameter represented a 
current density of approximately 1 A/cm’. This resulted in a material removal rate of 
0.5+0.01 um/s (CHEN and GERBERICH, 1988). As indicated in Fig. 5, the removal rate 
was very linear, although slightly higher at 0.6 m/s for a current density of about 
12 A/cm’. Based on this calibration, it is estimated that the depths, 1, from the 
fracture surface at which SACP were made are: ¢ = 1,2 and 5 wm +0.5 um. 

SACP were made at 25 keV with a spatial resolution of 5 um; it is thought that 
information in patterns obtained under this condition was being derived within 100 
nm of the surface (KACZOROWSKI and GERBERICH, 1986b). This depth is greater 
than that disturbed during the process of metal removal; therefore, the information 
obtained is essentially that from bulk material. The fact that the jet polishing surface 
area is small with respect to both the specimen surface and the size of the polishing 
bath dictates that temperature rise effects were small even at these high current 
densities. 

From the 36 locations in x, y and ¢, it was found that little difference could be 
found between the x distances greater than 200 wm from the outer surface. Thus, 27 
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Fic. 3. Fractography of a planar {100} fracture surface tested at AK = 13 MPa ./m under a frequency of 


3 Hz(1, 11), 0.3 Hz(II1) and 0.03 Hz(IV). Arrow indicates crack growth direction in this Fe-3 wt.%Si single 
crystal. 
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Fic. 6. Comparison of channeling line width changes in a calibration sample caused by fatigue deformation ; 
(a) shows the line widths at the locations shown by arrows in (b) and (c); (b) had 0 strain and (c) 
approximately 24% cumulative strain. 
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Fic. 7. Technique for measuring surface deformation involved direct measurement of the increase in grain 
boundary separation, 6, normal to the applied stress (arrows) : SACP were taken at (a) prior to deformation 
and (b) after deformation of the same grain. 
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a bn crystal {100}<010> of Fe-3wt.%Si loaded within the scanning electron microscope to 
; i _ Pa vm. (a) Secondary electron image of the crack tip region showing strong slip lines (dark 
- is contamination from liquid leaking from the crack tip) ; (b) Mohr’s circles of strain at low resolution, 
covering most of the region shown in (a); (c) Mohr’s circles of strain within a small region very near the 
crack tip. 
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Fic. 14. Selected area channeling patterns (SACP) at (a) 2 um, (b) 10 wm and (c) 200 um below the surface 
of a growing crack ; (d), (e) and (f) are microdensitometer traces of the same type of line in each of (a), 
(b) and (c). 
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of the 36 measurements were thought to be representative of plane strain, and the 
other 9 of plane stress. A slice of material 800 um thick was removed from the side 
of the specimen with a wire cutter in a spark machine, as shown in Fig. 4. This surface 
was then electropolished and SACP were made at 8 locations below the fracture 
surface of 10, 20, 50, 100, 200, 500, 1000 and 2000 um. The error in location of the 
two SACP nearest the fracture surface were + 1 ym, while the uncertainty in location 
of those farther from the fracture surface was within + 10% of the distance. 

Loading frequency was found to have little effect on SACP appearance. This was 
consistent with the previous finding (KACZOROWSKI ef al., 1986) that the total 
dislocation density at these same depths below the crack plane, as determined by 
transmission electron microscopy, was the same, independent of crack growth 
frequency. Thus, all the SACP taken at the same depth below the fracture surface 
were considered together. 


2.4. Analysis of channeling patterns 


Changes in SACP caused by varying dislocation density were quantified by meas- 
uring the width of specific lines on each pattern using an optical densitometer. This 
technique may be used to yield consistent results because the SACP were made from 
single crystals, thereby preventing problems with using different lines on different 
patterns, as may pertain to polycrystalline samples. The channeling line and location 
along that line must remain constant for this measurement technique to yield useful 
results. Figure 6 shows a typical SACP and the line and location used for line width 
measurement. SACP line width variations were converted into plastic strain by the 


use of a calibration standard. Previously, a unidirectionally loaded single crystal 
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Fic. 5. Penetration depth achieved in jet polishing an approximate 1.5 mm circular area for various times. 
Current density is ~ 1.2 A/cm’. 
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calibration sample had been prepared (KACZOROWSKI et a/., 1985), and the widths of 
the same line used for measuring differences from the fracture sample were correlated 
with strain, as measured from an extensometer mounted on this compression sample. 
Also used in this calibration were SACP of about the same orientation from poly- 
crystalline samples of a ferritic steel having about the same yield strength as the Fe- 
Si single crystal (KACZOROWSKI and GERBERICH, 1986). Comparable but slightly more 
consistent results were obtained in the polycrystalline samples, partly due to more 
uniform deformation and partly due to the accuracy of measuring local plastic strain. 
That is, the local plastic strain could be followed by measuring the separation of grain 
boundaries. Even though there was some nonuniform deformation from grain to 
grain, the same grain was followed from strain level to strain level by measuring the 
change in grain boundary separation, an example of which is shown in Fig. 7. Since 
this strain calibration was considered to be more accurate and an identical calibration 
on fatigue samples had been made on the same ferritic steel, this was considered to 
be the best available calibration curve. 

The calibration sample used for fatigue was obtained from a coarse grained, ferritic 
steel having about the same yield strength as the Fe—Si single crystal. The load level 
used was sufficient to produce plastic strain, and the number of cycles used at each 
strain increment to achieve a certain accumulated strain ranged from 50 to 330 cycles. 
SACP were made for each increment of strain and line widths were measured. The 
strains used for the calibration were measured in an identical manner to those for the 
monotonic calibration as indicated in Fig. 7. This procedure was believed to simulate 
the strain history reasonably well experienced by a material point just ahead of the 
crack tip as the crack grew toward it. For example, a material element ahead of the 
crack having a plastic zone size of 250 4m ahead of the tip would experience about 
125 cycles if the crack was growing at 2 um per cycle, see Fig. 8. From the monotonic 
and fatigue calibrations, the normalized line widths, w;/w,, were plotted as a function 
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Fic. 8. Schematic illustrating how points | and 2 move in the zone of plasticity as the crack grows. Points 

were initially outside the plastic zone boundary at (a); Point | moved to the tip with crack growth, while 

Point 2 moved through the plastic zone as the crack passed. The number of cycles over which strain 
accumulated for these two points is approximately given by R,/(da/dN). 
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of plastic strain. Here, the original line width, w,, refers to the measured value prior 
to loading. In Fig. 9 it is seen that line width increases with deformation as expected 
(CoaTEs, 1967). However, the amount of deformation to achieve the same line 
width increase is greater for the fatigue situation due to dislocation annihilation and 
formation of low energy configurations. This is discussed in greater detail elsewhere 
(KACZOROWSKI and GERBERICH, 1986). 

Having this calibration, it was possible to analyse strains from the single crystal 
crack propagation specimens by comparing SACP made at various locations to the 
calibration patterns to obtain w;/w,. Two sets of calculations were made from the line 
width observations. The first was associated with a line width minimum determined at 
the start of the study to be 2.55 mrad. The second set used an average valye since it 
was found that the minimum line width for the same line type was somewhat lower. 
An average of five readings gave a base-line to be 2.40+0.16 mrad. Since the cause 
of this change in minimum line width is unknown, both sets of data are considered 
valid and are reported. 


2.5. Stereoimaging 


One specimen was prepared specifically for detailed crack tip analysis by using the 
stereoimaging method to measure displacements (WILLIAMS et al., 1980). A crack was 
grown in this specimen on the (100) in the [010] in hydrogen by cyclic loading up to 
AK = 12.2 MPa ./m, R = 0, at 3 Hz. The specimen was then transferred to a special 
cyclic loading stage which fit within the scanning electron microscope (DAVIDSON and 


Naay, 1978). Loads were applied in steps between K,,;, = 5 MPa ./m and K,,,.x. The 
crack was observed to open progressively towards the tip as the load was increased, 
with an opening K;, of 13.0 MPa ./m. Finally, the load was increased to the maximum 
available, K, = 25.2 MPa, ./m, in an effort to obtain unstable crack growth. 


m= 1.05+0.059-€p 


m= 1.0 + 0.0836: €°7 
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Fic. 9. Calibration curves relating the relative increase of electron channeling line width (w;/w,) to the 
amount of tension or fatigue deformation. 





TABLE 2. Analytical strain distributions for stationary and growing cracks 
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Plane strain 
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3. ANALYTICAL AND NUMERICAL STRAIN DISTRIBUTIONS 


The results of analytical and numerical studies seeking to predict the strain dis- 
tribution within the plastic zone are given in Table 2. The formulae given in the table 
were obtained from the analytical and numerical references cited in Table 1. These 
predictions will be compared with experimental results in the next section. 

Since all predictions of strain distribution are in terms of the plastic zone size, R,, 
it is first necessary to determine this quantity. Equations (4) and (5) were used to 
compute R, for both stationary and growing cracks, since it has been shown that the 
coefficients used in these equations are approximately the same regardless of crack 
wake effects (DRUGAN ef al., 1982). All formulations in the table are in terms of 
plastic strain, ¢,. This is appropriate to the data obtained from analysis of SACP, 
since changes in the pattern line widths are caused only by dislocations. Conversely, the 
strains determined by stereoimaging are elastic plus plastic ; therefore, it is necessary to 
add +1 inside the brackets in the formulae shown in Table 2 for the comparison to 
be on an equal basis. 

For stationary cracks, the equations for plane stress and plane strain which include 
strain hardening, Table 2, Eqs (8) and (9), are the standard Hutchinson, Rice and 
Rosengren (HRR) relations. For growing cracks without strain hardening, the analy- 
sis shown in Table 2 for plane stress is a mode I analogy of the mode III solution 
(McCLINTOCK, 1958). For plane strain, the solution used is by RICE et al. (1980). 

For growing cracks with strain hardening, only the analyses of GAO and HWANG 
(1981) are available, Eqs (12) and (13). Their formulation includes the term 
[In (A/r)]"'~™. In order to be consistent with other formulations, we have assumed 


A = R,. Only the constants « and f are undetermined in Table 2. For the stationary 
crack, « may be determined from the HRR formulation with N = 1/n = 2.63 for 
power-law hardening, and for plane strain by using 


8, = by kr NN+ De} (14) 


2/(N+ 1) l—y? 1/(N+ 1) 
lal 


From a tabulation for plane strain results, /(2.63) is interpolated to be about 5.6. 
Working through (14), it can be shown that a, ~ 1/2 for 6 = 2/2. Since for a number 
of 0 values appropriate to this study, 1/2 < « < 1 for the stationary cracks, we use 
a; ~ | for all cases including «3. Finally, the f coefficient for the crack growing under 
plane strain has been determined to range from about 4.28 (DEAN and HUTCHINSON, 
1980) to 5.81 (LAM and MCMEEKING, 1984) with the most accepted value 5.46 
(DRUGAN et al., 1982). This latter value was used for both f and f’ in Table 2. It is 
realized that the «s and fs may be somewhat different than what were chosen. 
Although this might shift a given strain distribution up or down the strain axis, it will 
not change the character of the distribution. Given these caveats, a comparison of 
these models is made to the experimental strain distributions. 
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4. RESULTS AND DiSCUSSION 


4.1. Stationary cracks 


Stereoimaging was used to measure displacements in the near field of the crack tip 
for two conditions : the load required to just open the crack to the tip, K; = 13.0 MPa 
./m, and at the maximum applied load, K, = 25.2 MPa ./m. The corresponding loads 
were p = 1252 and 2424 N. The displacements caused by loading to the maximum 
value are shown in Fig. 10. The displacements, referenced to a point ahead of the 
crack tip, are relatively symmetric about the crack line. Strains derived from these 
displacements are shown in Fig. 11, where two relatively symmetrical shear bands are 
seen emanating from the crack tip. These results are consistent with the slip lines 
found on the specimen surface as shown in Fig. 12. Here, the crack tip with a large 
amount of surface plasticity is shown in Fig. 12(a). Mohr’s circles of strain for the 
overall crack tip region shown in the photograph are shown in Fig. 12(b) while Fig. 
12(c) shows the near crack tip region with higher resolution. In these figures, the 
diameter of the Mohr’s circle represents the maximum shear strain. Note that the 
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Fic. 10. Displacements around the crack tip shown in Fig. 12 at K, = 25.2 MPa /m as measured by 

stereoimaging. The crack tip is located at the cross (+), and loading was vertical. Note that the displacement 

calibration (right scale) is different from the spatial scale (left). Displacements behind the crack tip are the 
result of both crack opening and material deformation. 
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Fic. 11. Distribution of maximum in-plane shear strains around the crack shown in Fig. 12. Strains were 

derived from the displacements shown in Fig. 10. x and y are in micrometers ; loading axis was along the 

x-axis ; the crack is shown schematically on the plane of zero strain. Strain maxima may be seen along two 
bands at angles of 45 to 60° relative to the direction of crack growth. 


spatial and strain scales (left) differ from those at the right. Near-tip strain asymmetry 
was more accentuated by measurements made at higher resolution, Fig. 12(c). As 
there appeared to be a strain maximum close to a 60° ray, this was used to compare 
the theoretical models. To smooth out the strain asymmetry, values symmetrically 
positioned across the plane of the crack were averaged with “‘error-bars”’ representing 
the extremes. A similar procedure was followed for measurements made at the lower 
load. 

Tabulated data of &,.., &,, Yx» ANd Ymax, available at each point in the near crack tip 
region, were used to compare to the theoretical and numerical models using the above 
averaging procedure. Values of y,,,, versus distance from the crack along a 60° ray 
are shown in Fig. 13. From Eq. (8) in Table 2, and the calculated plastic zone sizes 
of 614 and 8500 um from (5) and (4),, the distribution of total strain was calculated 


from 
aes Oy R, 0.725 
‘7; | : | , (16) 


which is the total strain equivalent of (8). Here, the exponent comes from I/(n+ 1), 
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Fic. 13. Distribution of maximum shear strain at +60° to the crack line for the crack of Fig. 12 loaded 

to K = 13 and 25.2 MPa ,/m. Vertical bars represent the variation in strain along the two lines ; the circle 

marks the average value. The lines are from Eq. (8) of Table 2 using plastic zone sizes of 614 and 8500 um 
respectively. 


where, from (3), n = 0.38. The strain distributions computed using (16) are shown 
compared to experimental data in Fig. 13, which omits the strain at the crack tip. It 
should be noted that the crack tip strain is finite, while (16) predicts that it is not. 


Conversely, the strain distribution within the plastic zone is fit fairly well by the 
experimental data to a point within 2 ym of the crack tip. A least squares fit of these 
two sets of data gave slopes of —0.645+0.006 not too different from the theoretical 


slope of —0.725. 


4.2. Growing cracks 


As outlined in the experimental procedures section, distributions of strain were 
derived using SACP from both the surface and interior regions of specimens 
which had been broken by fatigue crack growth. Three SACP taken at 2, 10 and 
200 um below the crack plane are shown in Fig. 14, together with the correspond- 
ing microdensitometer traces. These data together with 91 other SACP were then 
converted to strain using the calibration curve in Fig. 9. 

Strain data so derived from SACP made at locations x, and x, at various y; in Fig. 
4 for a fatigue crack growing at AK = 13 MPa ./m are shown in Fig. 15. The 
horizontal bars represent uncertainty in location of the SACP, as described in Section 
2, while the vertical bars represent the spread in triplicate readings of the line traces 
and the uncertainty in the base-line calibrations. For K,,,, = 14.4 MPa ./m, the plastic 
zone size calculated using (4), of Table 2 was 750 um. This would appear to be a good 
estimate since no plasticity was encountered at a distance greater than 1000 wm. The 
solid lines in Fig. 15 were calculated using this plastic zone size in analyses for a 
growing crack, Eqs (8), (10) and (12) of Table 2. As may be seen from the figure, 
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Fic. 15. Plastic strain distribution perpendicular to the direction of crack growth for a fatigue crack grown 

at K,,., = 14.4 MPa ,/m in the near surface region (plane stress) as derived by SACP. Symbol © are based 

on a zero strain line width of 2.55 mrad, while squares are based on 2.40 mrad. Vertical bars are variations 

between locations sampled, and horizontal bars represent uncertainty in spatial location. Solid curves are 
Eqs (6), (8) and (10) of Table 2. 


none of these analyses fits the data well, although the elastic—ideally plastic case for 
a stationary crack gives a reasonable quantitative fit to much of the data. With strain 
hardening, the fit is improved, but still inadequate. The strain distribution for an 
elastic-plastic growing crack gives the best fit to the data, but underpredicts strain 
magnitude by a factor of three. 

The correlation of the three analogous models, (7), (9) and (11), to the plane strain 
distribution was also attempted. Here (5) gave a calculated value of 250 um for the 
plane strain plastic zone size. This also seemed reasonable since no strain was observed 
beyond 200 yum. As seen in Fig. 16, the correlation is not much better with either the 
shape of the curves being inadequate as in (7) and (11) or the strains being consistently 
low as in (9) for the stationary e—p” crack. For the latter, it is known from the HRR 
fields that «, < 1 and thus a better theoretical prediction would produce a greater 
discrepancy. It could be argued that the experimental error bars are large and do 
come close to traversing nearly all three of the models except at the largest strains. 
Pursuing this line of reasoning would indicate that such data are tenuous at best for 
testing such models. This plausibility is rejected for two reasons. First, if the maximum 
of the strain ranges is taken at each point, a curve through these gives a singularity 
not too different from the average and neither agrees with any of the theoretical 
models. Secondly, on the average there are about five measurements per point and 
the strain distribution in these single crystals tends to be localized in slip bands 
rather than be homogeneous. Even though reasonably homogeneous deformation was 
expected based upon the multiple slip of a large number of systems, on a scale 
commensurate with selected area channeling, slip was inhomogeneous. Thus most of 
the “error bar” is due to inhomogeneous deformation rather than experimental 
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Fic. 16. Plastic strain distribution perpendicular to the direction of crack growth for a fatigue crack grown 
at Kya, = 14.4 MPa ,/m in the specimen interior (plane strain) as derived by SACP. Symbol © are based 
on a zero strain line width of 2.55 mrad, while squares are based on 2.40 mrad. Vertical bars are variations 
between locations sampled, and horizontal bars represent uncertainty in spatial location. Solid curves are 
Eqs (7), (9) and (11) of Table 2. 


technique. It is the average then that should be representative of a continuum descrip- 


tion. 
Consider then the final correlation, a growing crack with strain hardening, which 


was most appropriate but most uncertain (HUTCHINSON, 1981, 1987). As stated above, 


the same pf’ = f = 5.46, R, = 251 um, n = 0.38 and (13) was used for plane strain. 
For plane stress, the parameters were «; = 1, R, = 753 um and n = 0.38 with (12). It 
was surprising to see how good the fit was in Fig. 17. Here, the “‘error bars” are the 
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Fic. 17. Data of Figs 15 and 16 compared to the theoretical plane stress and plane strain solutions of Eqs 
(12) and (13) for a growing crack with a work-hardening material. Only the most probable strain value at 
each location has been used for this comparison. 
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same as in Figs 15 and 16 but have been left off for clarity. Outside of the proper 
ordering of plane stress and plane strain, there is also a good quantitative fit at 
distances | um and greater away from the crack tip. If this or similar models result 
in additional corroboration, this implies that predictive models of ductile fracture in 
strain hardening materials would be on a firm basis. 

It should be mentioned here that after this correlation had been found, Fig. | 
was revisited. Here, the previously published parameters, namely, o,, = 218 MPa, 
E = 2x 10° MPa, n = 0.31 and R, = 1000 um were used with (12) from Table 2. This 
gave the dotted curve in Fig. 1. Note that the second form for (1) is similar to the 
calculated logarithmic singularities predicted for a growing crack except for a cut-off 
parameter, B, at r = 0. Here, B was taken to be one micrometer. The strains very 
close to the crack tip were relatively insensitive to this cut-off parameter. Here, 
computed strains would only vary by a factor of two at the crack tip for a cut-off 
parameter ranging from | nm to | wm ina sample with a 1000 um plastic zone. 


4.3. Other issues 


During the course of this investigation several important issues other than the 
distribution of strain became apparent. The first of these was that a growing crack 
appeared to have a larger accumulated strain magnitude than the static crack at a 
given position away from the crack tip. This might not be too unexpected since strain 
is accumulating as a crack grows under a material point. Comparing the plane stress 
results of Fig. 15 for the growing crack to Fig. 13 for the stationary crack, the strains 
for the growing crack appear to be larger at distances ranging from | to 750 um. If 
the strain-hardening formulation is right, and one compares (12) to (8), for surface 
strains this would predict greater strains for a growing crack over the range 
5000 > R,/r > 1 which corresponds to distances of 0.15 to 750 um. A similar cal- 
culation for plane strain using (13) and (9) gives 1040 > R,/r > 0.48 or 0.24 to 
520 um as the range. Note that this calculation was accomplished for total strain so 
that unity had been added inside the brackets of the equations. 

A second issue concerns failure criteria. This is particularly interesting since Fe- 
3%Si single crystals are considered to be semi-brittle at room temperature. Strains at 
| wm from the crack tip (Fig. 15) of 0.48 were obtained for a growing crack under 
plane stress and values of 0.26 (Fig. 16) under plane strain conditions. This was at 
Ky"** = 14.4 MPa ./m. Crack-tip shear strains of 1.48 (Fig. 12c) were obtained in 
plane stress with no failure and corresponding normal strains were 1.4. These were 
obtained for the stationary crack loaded to 25.2 MPa ./m. If this strain is put into 
the constitutive relation 


o = ke" = 1140 MPa, 


this is greater than the macroscopic cleavage fracture stress for Fe-3%Si which had 
been determined to be about 630 MPa (GERBERICH et a/., 1984). This implies that the 
fracture toughness of a highly dislocated Fe-3%Si single crystal is greater than its 
annealed counterpart. It is also interesting that this flow stress is about as great as the 
plastically-constrained flow stress obtained at the lower K;, values in the plane strain 
region. Here with n = 0.38 > 0.2, LAM and MCMEEKING’s (1984) result is used for a 
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growing crack with strain hardening. They provided a numerical solution at 


r= 1.133 x 10° ?(K;/o,,)’. 


which is 27 xm at Kj"** = 14.4 MPa ,/m. This distance is about 10 times larger than 
the crack-tip opening displacement and about 10 times smaller than the plastic zone 
size. Using their results for o/o,, along 0=0°, we find o,, = Gm = 3.750,, = 
1110 MPa. Thus, the constrained flow stress is only slightly higher in the plane strain 
center at low K, than in the highly work hardened surface at higher K,. The impli- 
cation is that the normal stress is always on the verge of promoting cleavage as the 
crack is growing with K; > 13 MPa ,/m. 

A recent paper on crack-tip stresses for anisotropic, elastic-ideally plastic material 
provides for even greater deviation in computed and measured cleavage stresses (RICE, 
1987). This shows that the stresses developed may be even greater for the BCC single 
crystal compared to its polycrystalline counterpart. Fora {010} <101) oriented crystal, 
the vertical (o,,) and parallel (@,,) stresses were shown to be 17.9 and 8.9% greater 
than the Prandtl field on the average. Although we cannot say how this would translate 
when the calculation is done for a {100} <010) crystal, the suspicion is that the stresses 
would also be higher there. How such large stresses can be maintained without 
cleavage has been addressed from a dislocation standpoint in one way by AsHBy and 
EmBury (1985) although the dislocation shielding aspects (Rice and THOMSON, 1974) 
alluded to in the introduction cannot be discounted. 


5. CONCLUSIONS 


Strain distributions to within one micron of a stationary or growing crack can be 
determined for both states of plane stress and plane strain. Both electron channeling 
and stereoimaging provide sufficient accuracy to give good measures of the distri- 
bution. In Fe-3wt.%Si single crystals with a {100}<010) crack, the strain distribution 
agreed very well with the HRR small scale yielding solutions for a stationary crack. 
For a growing crack, the distribution did not agree with any known solution except 
one by GAo and HWANG (1981) which includes strain hardening. This requires 
corroboration. Other issues such as the large magnitude of the stresses compared to 
known cleavage stresses in such brittle single crystals need to be addressed in terms 
of material resistance or dislocation shielding. 
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ABSTRACT 


RESIDUAL stresses in heterogeneous materials may arise because of differential or anisotropic thermal 
expansion of constituents. The paper is concerned with thermoelastic solids whose material properties 
fluctuate on the microscopic scale. Rigorous general relations between stored elastic energy and statistical 
averages (mean values and fluctuations) of residual stresses are derived. These results are applied to two- 
phase composites and to materials where the fluctuations of elastic constants can be neglected. One obtains 
exactly the stored energy, certain conditional mean values and the covariance matrix of the residual 
stresses. Under the assumptions of statistical homogeneity and isotropy, the results hold for any type of 
heterogeneous microstructure. 


INTRODUCTION 


WITHIN THE framework of linear elasticity we consider stresses in heterogeneous solids 
as polycrystals or composites. Due to differential or anisotropic thermal expansion 
of the components, residual stresses arise, which may influence mechanical properties 
of the material such as toughness or strength. It is our aim to derive some exact 
theoretical relations that can be applied to characterizing the fluctuating micro- 
structural stresses. This problem is closely connected with the calculation of the stored 
elastic energy, which deserves interest also from the thermodynamical point of view. 
For instance it may play a certain role in the thermodynamic free energy balance 
governing phase transformations. 

Usually the microstructure of heterogeneous materials shows a statistical disorder. 
Therefore the internal stress and strain fields fluctuate randomly, and it seems very 
difficult to derive general solutions for these fields. Accordingly, in the literature one 
can find theories which involve particular assumptions or simplifications as, for 
instance, two-dimensional models (Fu and Evans, 1985), special correlation functions 
of material constants (ORTIZ and MOLINARI, 1988), or special approximation schemes 
such as the effective medium approach (BOBETH and DIENER, 1987). Nevertheless 
there are several cases for which exact statistical expressions can be obtained, and it 
is the aim of the present work to find such results. 

The paper is organized as follows. After studying the basic equations in Section 2, 
we first derive some general relations which connect averages and fluctuations of the 
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random residual stress field with certain partial derivatives of the stored energy 
function (Section 3). The usefulness of these results is demonstrated by applying them 
to composites made of isotropic constituents (Section 4) and to materials where the 
fluctuation of elastic constants can be neglected (Section 5). Particularly in the latter 
case one obtains simple explicit results for the fluctuating residual stress and the stored 
energy that hold for every microstructural arrangement of components under the only 
assumption of statistical homogeneity and isotropy. 


2. Basic EQUATIONS 


We start by considering the relevant material constants. The linear thermoelastic 
behaviour is described by the constitutive law 


o=C(e—e') or ¢=So+e’,t (2.1) 


where o and ¢ denote stress and strain, respectively, C is the fourth-order tensor of 
elastic moduli, S = C~' is called the compliance tensor, and ¢’ is the tensor of stress- 
free strains (also called mismatch strains) which may especially arise because of 
thermal expansion. However other changes of shape or volume (e.g. due to phase 
transformations) may be comprised by e’. Thus 


» 
fix = | oi (T’) dT’ + &ix, (2.2) 
T 


0 


where « is the tensor of thermal expansion coefficients, 7) and T denote reference and 
actual temperature, respectively, and e" denotes a transformation strain. In this work 
the temperature is considered to be uniform over the solid whereas « and e" may vary 
according to the heterogeneous microstructure. So we shall regard ¢’ as the governing 
material property irrespective of its particular origin. 

In the heterogeneous medium, C = C(x) and e’ = e’(x) are certain random func- 
tions of position x. Consequently, the internal stress and strain fields become random 
functions, too. The appropriate description of random field variables is provided 
by expectation values such as mean values <o,(x)>, or second-order moments 
{0x (X)O,,(x)>, Where <...> means averaging in the statistical sense over an “‘en- 
semble” of samples. (Note: if an ergodic hypothesis applies, one can replace 
statistical averages by volume averages ; but this is not required in the present work.) 
As usual we restrict ourselves to statistically (i.e. macroscopically) homogeneous 
materials. Then an expectation value, <A(x)> say, of any field variable A(x) does 
not depend on position. 

Let us formulate the governing field equations. Equilibrium and compatibility are 
ensured by 


+ Bold face lower case greek symbols denote symmetric second-order tensors and bold face upper case 
letters are fourth-order tensors. We shall frequently use the common symbolic notation for tensor products : 


OF = OnE, Ce = CicimEims SC = SixtimCimpq 
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0 1/0 0 

~—Ox(X) =O, &(x) = ~~ U(X) + =— U(X) }, (2.3) 
Ox; 2 \ 0x; OX, 

where u,(x) is the displacement field. It is the aim of the present work to study the 
pure residual stress problem, hence we assume traction-free boundaries, 


O,.(xX)n(X) az () xeS, (2.4) 


where n, is the outward normal on the surface S of the body (superimposed external 
loading has been considered elsewhere, e.g. KREHER, 1988). 

Next we give some equations involving averages of stress and strain which immedi- 
ately follow from the basic equations (2.1, 2.3, 2.4). According to (2.4) the volume 
average of the stress vanishes, and, because of statistical homogeneity, this applies 
also to the statistical average : 


<a(x)> = 0. (2.5) 
The average strain is just the effective (i.e. macroscopic) stress-free strain e* : 
<e(x)> = e*. (2.6) 


In the special case of linear thermal expansion, s* is simply equal to the effective 
thermal expansion coefficients times the relevant temperature difference. 

The strain energy stored in the random residual stress field will play an important 
role in this work. Let us denote the average value of stored elastic energy per unit 
volume by u*. Then we have the equation 


<a(x)S(x)o(x)> = 2u*. (2.7) 


Furthermore we shall utilize another scalar equation namely the so-called Hill 
condition (HILL, 1963): 


(a(x)e(x)> = (a(x)><e(x)>.T (2.8) 
By (2.8), Eq. (2.7) with (2.1, 2.5) yields 
<a(x)e’(x)> = —2u*. (2.9) 


We note that the stored energy u* depends only on the fields C(x) and e’ (x) and is 
therefore a characteristic of the material. Hence u*, like e*, may be regarded as an 
effective material constant. This becomes obvious also by the fact that there are some 
close relations between u* and e*. Furthermore we remark that, in principle, u* can 
be determined experimentally by measuring the specific heat of the heterogeneous 
material (see e.g. CHRISTENSEN, 1979). 

For the considerations in the next sections it is sufficient to take into account only 
Eqs (2.5—2.9). These equations provide the relevant information determining mean 
values and fluctuations of the residual stresses. 


+ HILL (1963) proved (2.8) referring to volume averages. Due to statistical homogeneity, however, the 
Hill condition can be proved also for statistical averages (see e.g. KREHER, 1988). 
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3. GENERAL RESULTS 
It is our present aim to derive relations between the stored energy u* and certain 


statistical averages of the random stress. This will be achieved by means ofa variational 


procedure. 
We start by considering a small random variation of the elastic compliances while 


e’ (x) is kept constant: 
S(x) > S(x) + dS(x). (3.1) 
This implies corresponding variations of the internal fields : 
e(x) > e(x)+0e(x), o(x) > o(x)+060(x). (3.2) 
From (2.5) we obtain 
<da(x)> = 0. (3.3) 


Furthermore, de(x) and do(x) fulfil (2.3) so that the Hill condition (2.8) applies. Thus, 
with (2.5) and (3.3), we have 


(dae) = (da><e> = 0, (3.4) 
<ade> = <a><de> = 0. (3.5) 


(For simplicity the argument ‘(x)’ will frequently be dropped.) Let du*|,7 denote the 
variation of the stored energy caused by 6S(x) while e’(x) is kept fixed. Then, from 
(2.7) with (2.1), we find 


2du*|,7 = <adSa> +2<d0Sa> = <adSa) +2< dae) —2<da8'. (3.6) 


We note that (2.9) yields 
<da8’> = —26u*|,r. (3.7) 
Thus, from (3.6) with (3.4) and (3.7), one obtains finally 
<a0Sa> = —20du*|,r. (3.8) 
Next we consider a variation 
e’ (x) > e’ (x) +68’ (x), (3.9) 


while S(x) is kept constant now. Equations (3.2—3.5) apply analogously. We take the 
variation of (2.9) and find 


—2d6u*|, = <ade'>+<da8'> = <ade'> + (da8> — <(da(e—8')>. (3.10) 
By means of (2.1) and since dC = 0, the third term can be rearranged, 
da(e—s') = (e—8')do0 = (e—8' )Cd(e—8') = od(e—8') = ode—o0d8'. (3.11) 


Thus (3.10) becomes 
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—2du*|, = 2<ade8' > + <da8> — ade), (3.12) 

and because of (3.4) and (3.5), the last two terms vanish, and we find 
<ade'> = —du*|s. (3.13) 


Equations (3.8) and (3.13) can be evaluated in order to derive particular mean values 
and second-order moments of the residual stresses. This will be done in the next two 
sections by taking special variations dS(x) and de’ (x). 

We note that the case of internal stresses caused by external loading can be treated 
analogously. A relation corresponding to (3.8) was derived by BOBETH and DIENER 
(1986), and the analogue of (3.13) has been given by KREHER and POMPE (1989). 


4. COMPOSITES MADE OF ISOTROPIC COMPONENTS 


We consider a composite made of n components. Here S and e’ can take on the 


n 


values S, and ¢/ with probabilities p, (>) p, = 1, p, is usually called the volume 
r= | 


fraction). Thus the random fields S(x) and e’(x) may be written in the form 
S(x) = } S,O,(x),  #7(x) = ¥ ¢/0,(x), (4.1) 
r= | r=1 


where randomness is restricted to the characteristic functions ©, which describe the 
microstructural topology: 


1 if x € component r 


O,(x) = (4.2) 


0 if x € component r. 


We note the following relation: 
(A(x)O,(x)> = p,<A(x)|P>, (4.3) 


where <A(x)|r> denotes the conditional average of a certain field variable A(x) which 
is obtained when A(x) is averaged over a particular component r only. 

Let us suppose special variations of compliances and stress-free strains where the 
material properties of just one component change by 6S, or de but the material 
constants of the other components and the topology remain constant. Hence 


dS(x) = 6S,0,(x), de’(x) = 687O,(x). (4.4) 


Thus, considering (4.3), the results of (3.8, 3.13) turn into 


2 
(Fix Fim|1> OCSikim) = — p, du*|,7, (4.5) 


r 


l 
<Ox|r> O(E%), = P, du*|s. (4.6) 


r 


These general formulae will now be specialized for isotropic components where 
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l l / | ss 
Sikim >= Oi, 0m + 90cm + OimOx — 300m b) Ei ‘aon €; dix. 4.7 
(Sikim) —_— 4G, Ws iVki 30x 01m) (Eq) k ( ) 
Here K, G and e’ denote bulk modulus, shear modulus and isotropic stress-free strain, 
respectively. We presume that the stored energy is given as a function of component 
properties, 1.e. 


u* = a ( a a. a a ae m4) (4.8) 


Since we may vary K, and G, independently, we first take the variation 6K,. Hence 


Ay,* 


Se = : y . 
5) OKO, 01m, ou*|,7 eal OK,, (4.9) 
“ OK, 


O(Sixim), — OK 


and (4.5) yields 


; 
2K; Ou* 


4.1 
P, OK, ( 0) 


<s*|r> = 


(s = 6;/3, 6% = 8dy4+5,). Analogously, by taking variations 6G, and de/, we find 
from (4.5) and (4.6), respectively : 


4G? du* 
p, OG,’ 


1 ou* 


3p, Oe; 


(SiS ll) = (4.11) 


<s|r> = (4.12) 
In order to determine <s,|r> it would be necessary to calculate du* for a small 
deviatoric stress-free strain de),. However, in the special case where we can asssume 
that neither component of the composite prefers a particular orientation (statistical 
isotropy), the average component stress must be also isotropic (no direction can be 
distinguished). Hence 


<salr> = 0. (4.13) 


The results (4.10—-4.13) show that we can calculate mean values and second-order 
moments of residual stresses in the components of a composite in a simple way 
provided that we know the stored energy u* as a function of component properties. 
Let us demonstrate the usefulness of the theory by considering a statistically (i.e. 
macroscopically) isotropic composite made of two isotropic components. In this 
case the stored energy depends only on the difference of the two stress-free strains. 
Furthermore, u* must be proportional to the square of this difference. Hence 


u* = c(e5 —e')’, (4.14) 
where the factor c depends on the two bulk moduli, shear moduli and volume fractions. 
Thus we find 

A7,* 
= = ~206l~e) @ —2-—, (4.15) 
de} 


and (4.12) yields 
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u* 2 u* 


<s|1> = r—;, <sl2> = — —?. (4.16) 


3p; ese) 3p» es —e} 
Due to statistical isotropy, the effective stress-free strain has the form 

eX = e* Ox. (4.17) 
Hence (2.6) with (2.1) and (4.16) results in 


e* = Ce) = pi <ell>+p2<el2> = cep 2 a. ) T_,T (4.18) 


with 
<e"> = pie} +pre}. (4.19) 


Thus we have obtained a direct relation between the stored energy and the effective 
stress-free strain. Equation (4.18) can be rearranged: 


K, K, ~ 
o « 2 tte T_ ol) 4.20 
ay 3 (e* —<e"))(e} —e7) (4.20) 
(The case K, = K, will be considered in Section 5.) We remark that the result (4.20) 
is in accordance with the expression for the specific heat given by ROSEN and HASHIN 
(1970) (the latter was derived by means of bounding theorems). In addition we note 
that there is also a relationship between e* and the effective bulk modulus K* (Levin, 
1967, see also ROSEN and HASHIN, 1970): 
oa fe p,/K,+p.2/K,—1/K* oF . 421 

e ce >+ 1/K, -1/K, (e; —e}). (4.21) 
Thus one may use any theoretical result for the effective bulk modulus or the effective 
stress-free strain to calculate u* = u*(K,, K>, G;, G>, e|, e}). From this stored energy 

the residual stress fluctuations are found by applying (4.10—4.13). 


5. HOMOGENEOUS ELASTIC CONSTANTS 


Especially interesting results are obtained if the fluctuations of elastic constants are 
negligible. So let us assume 


aa D ietins al a ie 
Sikim(X) — Sikim = 0x0 im + 4G (00m + Dim Oki — 30i01m)- (5.1) 


9K 
The stress-free strains e’ (x) fluctuate. For the sake of simplicity we presume a discrete 
probability distribution, i.e. e’(x) can have only discrete values e/, and we adopt the 
description given in (4.1-4.3) (later on we shall take into account continuous dis- 
tributions also). 

As an example one may imagine a one-phase polycrystal where the elastic anisotropy 
of single crystals can be neglected (this applies to Al,O,, for instance). In this case 
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es 


the anisotropic thermal expansion in combination with a discrete orientational dis- 


tribution of crystallites leads to a discrete set of values ¢/. 
In the present case we suppose that the stored energy is known as a function of the 
common elastic constants and of the individual stress-free strains: 


(5.2) 


The calculation of mean values and second-order moments of the random residual 
stresses can be accomplished just as in the preceding section. Instead of variations of 
individual compliances S,, however, it is now only possible to consider the variation 
of the common compliance S. Accordingly one finds just the total second-order 
moment of stress, and we have (cf. (4.5, 4.6)), 


(Cig Tim> OSikim _ — 206u*|,7, (5.3) 


scp Pia 
(ix TD O(E% ), — — du*|s. (5.4) 
Considering (5.1) and taking independent variations of K and G, (5.3) can be evaluated 
analogously with (4.9-4.11). The result is 
, Ou* , Ou* 


aK? SSuSa> = 4G? =. (5.5) 


(s*> = 2K? 
Furthermore we may take independent variations of the tensor components of ¢}, 
and we obtain from (5.4) (cf. 4.12): 
1 du* 
(ox \r> = — T° (5.6) 
P, O(Eix), 
It is worth noting that in the present case of homogeneous elastic constants, the 
effective stress-free strain is simply given by 


e* = <8’) = y p,st. (5.7) 


r= | 


This is easily seen from (2.5, 2.6), with (2.1). 
Now we evaluate (5.5, 5.6) in the special case of a statistically isotropic material. 
As shown in the Appendix the stored energy can be exactly calculated : 


— = 5 Ehim< (Ek —€&%>)(Eim — <Eim>)>- (5.8) 


The tensor E’ is given by (A14). Thus by decomposing ¢’ into its isotropic and 
deviatoric part we may alternatively write 


I8KG 9K+8G 
u* = 3K14G «© —<e’ >)?>+ tG 3K+ 4G <eieeik>? (5.9) 


(because of statistical isotropy we have <e},> = 0). Inserting (5.9) into (5.5) yields 
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) I2KG Y — 
(s*) = (en) (<(e7 —<e"))*) +40 <enen>), (5.10) 


2 ee _, 27K? +48KG+32G? | 
Jeter e146" GKsagy? eee = G.I) 


3 12KG 
(Sit Si > = 2\3K+4G 


and, from (5.6) with (5.8), we obtain 
(O%\r> = “im( (Ehime — <8im >) (5.12) 


Obviously the results (5.10—5.12) hold also for continuously varying stress-free strains 
(i.e. instead of the discrete probabilities p, we have a probability distribution p(e’)). 
Thus we may generalize (5.12): 


(Oi |e” —S ‘elm (Elm cr <Eim>)s (5. I 3) 
or, with (A14), 


12KG 


_ i ad jin =— 
3K14G <é >) ik . 


7 
(Ox |@ > = 
Here <o,,\e’ >) means the average value of residual stress for a particular region of the 
heterogeneous solid characterized by the stress-free strain e’. 
Due to statistical isotropy the tensor <o,¢,,,> must be an isotropic tensor. Hence 
it can be expressed by the two scalars <s?> and <s,5, >: 


(Onn Fim> - <s*> Oi Om ? io (SpqSpq > (OuO Km + Dim ki = 301 Sim): (5. I 5) 
So, from (5.10, 5.11), we find 
{Fim > = 4G? (Riki (e” — Ke")? + Rin Cpa pa >) (5.16) 


5 


meer a! |. Peer eee 
Rikim = 5 (BK 44G)? PuPkm + ims + Fix im), (5.17) 


ne. 27K* +48KG+32G’ “—s oe ve. K? 5.4 
ikim ~— 50 (3K+4G) (< il km im“ kl 3% ik“ Im 5 (3K+4G)? ik“ lm» 


(5.18) 
or, by introducing Poisson’s ratio v, 


| L/itvy.. — or 
Rikim _ 15 ( ad ) (OO Km + Oim OKI + 60 i, Om)s (5.19) 


1 19—34v4+19v* _ aos 1 1—6v+v? 
2? (00 Km + OimOxi) bal 75 ( | a vy)? 


Sixdm- (5.20 
450 (1—v)? we — 


. _— 
Riki — 
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It should be emphasized that the results (5.9, 5.14, 5.16) are valid for any het- 
erogeneous microstructure. Besides statistical homogeneity and isotropy, we have 
only assumed homogeneous elastic constants. The stress-free strain e? may randomly 
vary in an arbitrary manner. Therefore (5.16) provides exactly the total fluctuation 
(also called covariance matrix) of the random residual stresses (remember <a> = 0), 
and (5.14) gives the conditional mean value of the stress in a particular component 
or region characterized by the local stress-free strain e’. 

We remark that the validity of the results can be checked by introducing (5.14) and 
(5.16) into (2.9) and (2.7), respectively. One arrives again at the result (5.9) for the 
stored energy as it must be. 

The error introduced by the assumption of homogeneous elastic constants can be 
assessed by considering (5.14, 5.16). It turns out that the relative error in the pre- 
dictions for average stress and fluctuation of stress is of the order of the relative 
variation of elastic constants. Therefore the present results should be applicable if the 
latter variation is less than the required accuracy of stress prediction. Otherwise one 
might use the equations in the sense of upper or lower bounds on the stresses. 

For illustration it is useful to consider two special materials. First we evaluate the 
results for a composite made of two isotropic components (¢, = e76,). Here (5.9) 


turns into 


i I8KG te 
u = PiP23K 4 4G (2 -e) > (5.21) 


and (5.14, 5.16) yield 


I2KG 


P| 3K+4G (e5 —e7)ox, (5.22) 


rae 
<x |1> = pr 3K +4G (e3—e})b%, <o%|2> = — 


- | eae ea i aa tat Se 
<i Fim > -_ sapiva( Xe) (0/0 km + 0imOK1 + 665m) (C5 —e71)?. (5.23) 


Note that (5.21) represents the exact limit of (4.20) for the case K, = K,. The results 
(5.21 and 5.22) are in accordance with previous results (KREHER, 1988) which were 
derived by means of a direct limiting process. 

As a second example we consider a one-phase polycrystal. Here e’ varies only 
because of the different orientations of crystallites. The scalars e’ and e/e! do not 
depend on orientation, therefore we can drop the averaging brackets in (5.9, 5.14 and 


5.16). We obtain 


eee. . 


3K +4G 6% deh i | a el el (5.24) 


*—I1¢ 
u 5 pq™ pq? 


9K+8G , isa 


T 2 
i 2 = —<zG ci. 


= 2 pid) {7 pT 
(Fix Fim)> i 4G Rika Cpe pes (5.26) 
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with R“ given by (5.18) or (5.20). These results hold in any coordinate system, yet it 
is especially useful to adopt the crystal coordinate system where 


yy ee 
c=10 2 0 (5.27) 
0 O @¢ 
Then we have 
ener, = 4((e1 —05)? + (3 — 04)? + (e5 —e1)’). (5.28) 


The results (5.24 and 5.25) are in accordance with the special results for transverse 
isotropy («| = ¢5) given by KREHER (1988). 

As a practical example let us analyse the residual stresses in polycrystalline zirconia 
(ZrO,). SCHMAUDER et al. (1984) report single-crystal thermal expansion coefficients 
11) = %9 = 11.6 10° °/K, «3; = 16.8 x 10° °/K. In order to find the residual stresses 
arising during cooling from fabrication temperature down to room temperature, we 
assume a temperature difference 7— 7, = — 1500 K. Thus, from (2.2), the relevant 
stress-free strains are e| = ¢€3 = —17.4x10-°, ef = —25.2x 10°. We adopt iso- 
tropic elastic constants, K = 192 GPa and G = 88 GPa (this is the most usual assump- 
tion for this type of material). Then (5.25) yields 


in 0 0 
<ox\e'> = 2401 O —1 O| MPa. 
0 0 2 


This is the average residual stress of a certain crystal species, i.e. the mean value for 
crystals of a certain orientation (of course, due to statistical isotropy, the result does 
not depend on the particular orientation). In addition to this conditional average, the 
fluctuation of stress can be derived from (5.26). To get an impression about the 
magnitude of these fluctuations one may calculate scalar invariants. So, for instance, 
one finds 


SK Onon> = 630 MPa. 


These results can be applied to assess the risk of shrinkage microcracks, for example. 

Finally we compare our expression (5.26) with the corresponding result obtained 
by OrTIz and MO iNnarRI (1988). Under the assumption of a particular correlation 
function of stress-free strains those authors tried directly to calculate the covariance 
matrix (there denoted by <60;,60;,,>, but since <o,> = 0 one has do, = o,). The 
calculations are more tedious than the present ones and, moreover, must involve a 
mistake since the result predicts fluctuations which are about two times larger than 
those given by (5.26) with (5.18). This is directly seen by calculating the stored energy 
from (2.7) in combination with Ortiz and Molinari’s result for <o,¢;,,>. One finds 
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, ,IK+ 10G 7v 
uou = ; , Je Cf (5.29) 


Tt 
3K+4G rn it tet 


which does not agree with the exact value (5.24). In fact it is even larger than the 
general upper bound for u* which is given by (see e.g. KREHER and Pompe, 1959): 


ut = $¢e’Ce’) —h¢e'C)<(C) ‘Ce’. (5.30) 


In the present case of polycrystals made of elastically isotropic single-crystals, this 
upper bound reads 


‘ 7 7 
ul = Cf 4 < ur MM: (5.31) 


6. CONCLUSIONS 


The present paper shows that it is possible to derive exactly several expressions 
describing mean values and fluctuations of a random residual stress field. This can be 
done without particular assumptions about the microstructure. In the general case, 
the relevant microstructural information is provided implicitly by the stored energy 
which must be given as a function of component properties. This stored energy 
function may be obtained from other theories (in principle also by experiments). 
Particularly with two-phase composites, however, the stored energy can be derived 
from other macroscopic parameters such as the effective bulk modulus. 

Especially useful results have been found for materials where the fluctuations of 
elastic constants can be neglected. Here the conditional mean value and the covariance 
matrix of residual stresses and the stored energy are completely determined by the 
fluctuation of stress-free strains. These results do not depend on additional details of 
the heterogeneous microstructure as grain shape or special correlation functions of 
material properties. 


Note added in proof 
Ortiz and Molinari have corrected their formulae (see Corrigendum in this issue 
on p. 129). Now their results are in full agreement with the results in this paper. 
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APPENDIX: CALCULATION OF THE STORED ENERGY (5.8) 


At first we consider a heterogeneous solid with a given field of stress-free strains e’ (x) and 
homogeneous elastic constants C(x) = C. The basic field equations (2.3) together with the 
stress-strain relation (2.1) yield the following equation for the displacement field : 


CO? oO 
C iim « ~ u,,,(X) - C ikim « e},, (x) = 0. (Al) 
OX, CX; OX, 
Boundary conditions are prescribed by (2.4), i.e. we assume a traction-free surface. The solution 
of (A1) is given by 


( 
u(x) = -| 9G: (XX)C jem Ax’ Ein (X’) Ax’, (A2) 
x Ox 


ooR 


where B denotes the solid and g,,(x, x’) is the Green function of the medium with homogeneous 
elastic constants C and traction-free surface (see, for instance, DEDERICHS and ZELLER, 1973). 
By partial integration and taking into account the boundary conditions one obtains 


a 
Ey (Xx) = (; = uc) = | G iim (X, X') Cimpg®pg(X’) AX’ (A3) 
FT (ik) B 
with 
oC? 
G eim (X, x’) = C ‘ ~_/ G(X, x’) cimycim) + (A4) 
“Ak VAm 


The stress field is found by using (2.1) 
Oy(x) = -| Fim (XX )Eim (X’) AX’, (A5) 
B 


with the stress Green tensor defined by 
G ici (%; x’) -_ CictmO(X Tes x’) -* Ciepg Gpgst (Xs SC tus (A6) 


where 6(x) denotes Dirac’s 6 function. Due to the traction-free boundaries, the stress must 
vanish if we are concerned with a homogeneous field ¢’(x) = e”. Thus (A5) yields 
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| Goeim (X, X’) dx’ = 0. (A7) 
B 


Next we take into account the randomness of ¢7(x). Due to statistical homogeneity, <e’ (x)> 
does not depend on position so that (A5) with (A7) can be rewritten: 


Oi. (X) — -| tkim (X, x’) (Ein (X’) Sot (Ein >) dx’. (A8) 
B 


The stored energy u* may be calculated by means of (2.9). Due to (2.5) we can write 
u* = — Coq (x)EK(X)> = — 3G (X)(EK(X) — CER). (A9) 
By inserting (A8) we find 


| . 
u* = | G Sct (XX) K kim (X, x’) AX’, (A10) 
B 


with 
K ieim(XX’) = €(&(%) — (84>) (Eim(X’) — Cin >) >: (All) 


Due to statistical homogeneity, the correlation function K’ depends only on x—x’. We now 
adopt two additional assumptions about the random field e’ (x): 


(i) statistical isotropy, i.e. K’(x, x’) = K’(|x —x’}), 
(ii) no long-range correlation. 


The latter assumption is fulfilled in most heterogeneous solids. It means 
lim K’‘(x,x’) =0. (A12) 


Ix—x'|> ax 
Due to (A12) we can let the surface of the heterogeneous solid go to infinity, and we evaluate 
(A10) by using the infinite medium Green tensor. In this case, following KRGONER (1977), G” 
can be represented by two tensors E’ and F’, 


G’ (x, x’) = E’0(x—x’)+F* (A13) 


|x—x’|*’ 


bi A ok gi ee 
ae 3K 4G ik“ lm 5 3K+4G 


The tensor F’ depends on the direction of x —x’ and has the property 


(OO cm + SimOki = 30x Om): 


| ikim AQ = 0, (Al5) 


where dQ means an integration over all directions. 
Now we introduce (A13) into (A10). Since K’ depends only on |x —x’|, the term involving 
F’ vanishes and we find 


u* = DE ict K iim (Xs x). (A16) 


Thus, with (A11), one arrives at the result (5.8), which holds for any statistically homogeneous 
and isotropic material with homogeneous elastic constants but fluctuating stress-free strains 
> 

&° (x). 
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CORRIGENDUM 


MICROSTRUCTURAL THERMAL STRESSES IN 
CERAMIC MATERIALS 


M. Ortiz and A. MOLINARI 


J. Mech. Phys. Solids 36, 385-400 (1988) 


Equation (3.24) should read : 


— : . , 
te = 3 [(~a, —a2)° + (a, —%3)° + (42 — 3) *], 


9 . . , 
Be = 6 Pi —B)° +(B, —B3)° + (B2—B3)’]. 


Formulae (3.28) and (4.3) should be corrected to: 


2292 19v? —34v+ 19 7 gr 2v — 12v+2 
= gm > = Ae “6 5 
Ame aaa 9) 75(1—vy? 


32 yv*—v+1 


225 (1—v)? ° (4.3) 


g°(v) = 
The rest of (4.3) is correct. 

These corrections have an effect on the discussion on Al,O, given at the end of the 
paper. The computed value of the variance of the normal tractions changes to o = 
84 MPa, which is still in good agreement with the experimental range 80-100 MPa 
measured by BLENDELL and CosLe (1982). In addition, for grain sizes a = 50 and 
100 um, one computes fractions of microcracked facets f = 0:75 and 3.4% respec- 
tively. The corresponding effective moduli then become EF = 0.997E for a = 50 um, 
and E = 0.985E for a = 100 um. 

That the aforementioned formulae were indeed in error has been noted by KREHER 
(1990) in a recent paper. The above corrections bring our results into full agreement 
with Kreher’s. 


REFERENCES 


BLENDELL, J. E. and Coste, R. L. 1982 J. Am. Ceram. Soc. 65, 174. 
K REHER, W. 1990 J. Mech. Phys. Solids 38, 115. 


129 





J. Mech. Phys. Solids Vol. 38, No. 1, pp. 131-132, 1990. 
Pergamon Press plc. Printed in Great Britain. 


International Centre for Mechanical Sciences 


Centre International des Sciences Mecaniques 


I-33100 UDINE (Italy), Palazzo del Torso, Piazzo Garibaldi, 18 
Tel. (0432) 29 49 89, Fax (0432) 50 15 23 


PRELIMINARY PROGRAMME 1990 


Courses 
Modelling macroscopic phenomena at liquid boundaries 
Coordinators: W. Kosinski (Warsaw); 
A. I. Murdoch (Strathclyde) 25-29 June 
Reliability problems: general principles and application in mechanics 
of solids and structures 
Coordinators: F. Casciati (Pavia); J. B. Roberts (Brighton) 2-6 July 
Engineering applications of dynamics of chaos 
Coordinators: W. Szemplinska-Stupnicka (Warsaw); 
H. Troger (Wien) 9-13 July 
Shape and layout optimization of structural systems 
Coordinator: G. Rozvany (Essen) 16-20 July 
Linear prediction theory 
Coordinators: L. D. Davisson (Maryland); 
G. Longo (Trieste) 23-27 July 
Mechanics of robots 
Coordinator: J.C. Guinot (Paris) 3-7 September 
Stability problems of steel structures 
Coordinators: M. Ivany (Budapest); M. Skaloud (Prague) 24-28 September 
Progress in computational analysis of inelastic structures 
Coordinator: E. Stein (Hannover) 1-5 October 
Diagnostics of machinery 
Coordinators: C. Cempel (Poznan); O. Mahrenholtz (Hamburg) 8-12 October 


Other Events 
[IUTAM-CISM Course on ““Waves in Fluids” 

Coordinator: M. J. Lighthill (London) 17-22 May 
CISM’s 20th Anniversary 23-24 May 
8th CISM-IFToMM Symposium on Theory and Practice of 
Robots and Manipulators (Cracow) 

Chairmen: A. Morecki (Warsaw); G. Bianchi (Milan) 2-6 July 
CISM-IIASA Summer School on Methodology, Implementation 
and Applications of Decision Support systems 

Coordinators: A. Lewandowski (Wien); 

P. Serafini-M. G. Speranza (Udine) 17-21 September 
Consolidating Course on Fundamentals for Mechanics (Budapest) 
Coordinator: P. Scharle (Budapest) 1-9 October 





Announcement 


DESEG—Department of Structural Engineering and Geotechnics 
(in Italian ) 
Corrosione e protezione delle armature in c.a. 
Coordinator: P. Pedeferri (Milano) 
I] sistema di posizionamento globale satellitare (GPS): Metodologie 
ed applicazioni all’ingegneria del territorio 
Coordinators: F. Crosilla (Udine); L. Mussio (Torino) 


Gallerie 

Coordinators: K. Kovari (Zurich); M. Mele (Roma) 
La diagnostica ed il controllo di qualita delle opere di ingegneria 
civile con i metodi di indagine non distruttiva 

Coordinators: M. Mele (Roma); D. Almesberger (Trieste) 


Meetings Hosted by CISM 
Workshop in Microgravity Fluid Physics 
Coordinator: J. M. Haynes (Bristol) 


March 


3-5 April 


May 


7-9 November 











CONCISE ENCYCLOPEDIA OF 
WOOD & WOOD-BASED MATERIALS 








Edited b 

ARNO P SCHNIEWIND, 

Protessor of Forestry, Forest 
Products Laboratory, University 
of California, Berkeley, USA 





The Concise Encyclopedia of Wood & Wood- 
Based Materials is a volume in the new series 
Advances in Materials Science & Engineering, 
a new series of concise encyclopedias in 
major areas of materials science which 
combine newly commissioned articles with 
original material from the Encyclopedia of 
Materials Science & Engineering 


Contact the Marketing department at your 
nearest Pergamon office for more information 
on other titles in this series. 





a 


June 1989 
150 illus 
US$125.00 











Approx 372 pp 800 refs 
0 08 0347266 Hardcover 











Available in the Americas from MIT Press, 
Cambridge, Massachusetts, USA. 


Prices and publication dates are subject to change 
without notice. 


PERGAMON PRESS 


Member of Maxwell Macmillan Pergamon 
Publishing Corporation 





PERGAMON PRESS plc 

Headington Hill Hall 
Fairview Park Oxford 
Elmsford, NY 10523 OX3 OBW 
USA UK 


PERGAMON PRESS Inc. 
Maxwell House 


Alphabetically organized, with 73 articles by over 60 leading authorities on 
wood, the Concise Encyclopedia of Wood & Wood-based Materials covers 
the whole range of knowledge and current research in wood science, 
including articles on the availability and economics of timber resources, 
wood products such as plywood and mineral-bonded wood composites, 
descriptions of the major commercial wood species of the world, 
fundamentals of wood properties and behaviour, factors causing 
deterioration and their control, and principal processing methods. 


The Concise Encyclopedia of Wood & Wood-Based Materials draws relevant 
articles from the highly acclaimed Encyclopedia of Materials Science & 
Engineering. Where necessary, a complete revision of the original article has 
been undertaken. New articles have also been commissioned to provide 
up-to-date accounts of new developments in the area. 


Extensively illustrated, with nearly 200 photographs, drawings and tables in 
over 300 pages, each article is intended to serve as the first source of 
information on a given topic. The reader is guided to further reading by 
helpful cross-references and nearly 500 up-to-date citations in the 
bibliographies at the end of each article. A comprehensive, three-level 


subject index is also provided. 


The Concise Encyclopedia of Wood & Wood-based Materials will be 
invaluable to architects, engineers, builders, plant managers, wood 
technologists, purchasing agents, timber merchants. 





A SELECTION OF ARTICLES 











An Introduction to Wood and 
Wood-Based Materials 

Acoustic Emission and Acousto- 
Ultrasonic Characteristics 

Acoustic Properties 

Adhesives and Adhesion 

Archaeological Wood 

Bamboo 

Biotechnology in Wood Processing 

Building with Wood 

Cellulose: Chemistry and 

Technology 

Cellulose: Nature and Applications 

Chemical Composition 

Chemically Modified Wood 

Chemicals and Liquid Fuels from 
Wood 

Coconut Wood 

Constituents of Wood: Physical 
Nature and Structural Function 
Cork 

Decay During Use 

Deformation Under Load 

Density and Porosity 

Design with Wood 

Deterioration by Insects and Other 
Animals During Use 

Drying Processes 

Electrical Properties 

Energy Generation from Wood 

Fire and Wood 

Fluid Transport 

Forming and Bending 

Fumigation 

Future Availability 

Glued Joints 

Glued Laminated Timber 

Hardboard and Insulation Board 

Health Hazards in Wood Processing 

History of Timber Use 

Hygroscopicity and Water Sorption 

Joints with Mechanical Fastenings 


Laminated Veneer Lumber 
Lignin-Based Polymers 
Lumber: Behavior Under Load 
Lumber: Types and Grades 
Machining Processes 
Macroscopic Anatomy 
Mineral-Bonded Wood Composites 
Nondestructive Evaluation of Wood 
and Wood Products 
Paper and Paperboard 
Particleboard and Dry-Process 
Fiberboard 
Plywood 
Preservative-Treated Wood 
Protective Finishes and Coatings 
Radiation Effects 
Rattan and Cane 
Resources of Timber Worldwide 
Stringed Instruments: Wood 
Selection 
Shrinking and Swelling 
Strength 
Structural-Use Panels 
Structure, Stiffness and Strength 
Surface Chemistry 
Surface Properties 
Thermal Degradation 
Thermal Properties 
Timbers of Africa 
Timbers of Australia and 
New Zealand 
Timbers of Canada and the USA 
Timbers of Central and 
South America 
Timbers of Europe 
Timbers of Southeast Asia 
Timbers of the Far East 
Trade, Prices and Consumption of 
Timber in the USA 
Ultrastructure 
Weathering 
Wood-Polymer Composites 








EUROPEAN POLYMER JOURNAL 


Founder Editor: M STACEY, FRS, University of Birmingham, UK 

Joint Editors-in-Chief: JC BEVINGTON, Department of Chemistry, The University, Lancaster, UK and 
J V DAWKINS, Department of Chemistry, The University of Technology, Loughborough, Leicestershire, 
UK 


This journal acts as a medium for the exchange of research in the area of macromolecular substances, both 
synthetic and natural, and publishes results bearing on the physics or chemistry of polymers. Papers may 
fall either into the category of original papers dealing with substances of high molecular weight, or that of 
review articles covering the literature on topical aspects of polymer science and technology. The languages 
of the journal are English, French, German and Italian; for all papers an abstract is provided, and in the 
case of the latter three languages the abstract also appears in English translation. 


A Selection of Papers 

H DAOUST & P FERLAND (Canada), Solvent-induced conformational change of poly-(l-ornithine) and 
poly(l-lysine) hydrobromides in water-2-chloroethanol mixed solvent. 

R S DAVIDSON & C LOWE (UK), A study of some photocrosslinkable resins using i.r spectroscopy. 
P HODGE & A A NAIM (UK), Synthesis of crosslinked polymers containing benzoyl peroxide groups 
and their use as initiators. 

M NASIR & C H CHOO (Malaysia), Cure characteristics and mechanical properties of carbon black 
filled styrene-butadiene rubber and epoxidized natural rubber blends. 

A MATSUMOTO, H ANDO & M OIWA (Japan), Gelation in the copolymerization of methyl 
methacrylate with trimethylolpropane trimethacrylate. 

M A ZHURAVLEV & V B IVANOV (USSR), Mechanism of photochemical dehydrochlorination of 
poly(vinyl chloride). 


Indexed/Abstracted in: Curr Cont ASCA, Cam Sci Abstr, Chem Abstr Serv, Curr Cont/Phy Chem & Earth 
Sci, INSPEC Data, Pol Sci Tech Data, RAPRA Data, Curr Cont Sci Cit Ind, Curr Cont SCISEARCH Data. 


(00294) 

1990: Volume 26 (12 issues) 

Annual subscription (1990) DM 1420.00 
Two-year rate (1990/91) DM 2698.00 
ISSN: 0014-3057 


PERGAMON PRESS 


Member of the Maxwell Macmillan Pergamon Publishing Corporation 


Pergamon Press plc, Headington Hill Hall, Oxford OX3 OBW, UK 
Pergamon Press Inc., Maxwell House, Fairview Park, Elmsford NY 10523, USA 


Advertising rate card available on request. Back issues and current subscriptions are also available in microform. 
The DM prices shown include postage and insurance. For subscription rates in the Americas, Japan, UK and Eire please contact your nearest 
Pergamon office. Prices and proposed publication dates are subject to change without prior notice. 


























ENGINEERING FRACTURE MECHANICS 


An International Journal 


Editor-in-Chief: HAROLD LIEBOWITZ, School of Engineering and Applied Science, The George 
Washington University, Washington, DC 20052, USA 


Engineering Fracture Mechanics publishes technical papers on research and advanced applications of 
fracture mechanics, including contributions in the areas of mechanics and materials science as related to 
fracture mechanics. 

Topics regularly covered include: analysis of crack patterns, crack stress field analysis, experimental stress 
analysis methods, evaluations of materials, uses of fracture mechanics in studies of fatigue fracturing and 
stress corrosion cracking. Several pages of appropriate issues are devoted to a compendium of information, 
including a synthesis and utilization of significant research results. This information takes the form of 
reviews of books and papers published elsewhere, compilations of solutions (formulae) and other 
analytical results, and tables and curves of data and information from various sources. It provides a 
valuable data and reference source for the field. 


A Selection of Papers 

IS RAJU & W B FICHTER (USA), A finite-element alternating method for two dimensional mode I 
crack configurations. 

P S THEOCARIS (Greece), Definition of the crack tip as the center of the caustic. 

E SMITH (UK), The effect of a kink, a U-loop and a T-junction on the instability of a circumferential 
through-wall crack in a stainless steel piping system. 

X S ZHANG (PRC), Transient response of an edge crack in a rectangular plate with stress-free edges for 
the tearing mode. 

HUSEYIN SEHITOGLU & WEI SUN (USA), The significance of crack closure under high temperature 
fatigue crack growth with hold periods. 

D J QUESNEL & R W BUSH (USA), The destabilizing influence of plastic yielding on fracture. 

JS SHORT & D W HOEPPNER (USA), A global/local theory of fatigue crack propagation. 

Y NAKAJIMA, Y IINO & M SUZUKI (Japan), Fracture toughness behaviour of service-exposed type 
321 stainless steel at room and elevated temperature under normal and low straining rates. 


A Software Survey section is included in this journal. 


This journal contains a Cauzin Softstrip, providing each issue’s table of contents in a machine-readable 
form. 


Indexed/Abstracted in: Current Contents, Engng Ind. Monthly & Author Index, PASCAL-CNRS 
Database, INSPEC Database 


(00322) 

1990: Volumes 35-37 (18 issues) 

Annual subscription (1990) DM 2650.00 
Two-year rate (1990/91) DM 5035.00 
ISSN: 0013-7944 


PERGAMON PRESS 


Member of the Maxwell Macmillan Pergamon Publishing Corporation 
Pergamon Press plc, Headington Hill Hall, Oxford OX3 OBW, UK 
Pergamon Press Inc., Maxwell House, Fairview Park, Elmsford NY 10523, USA 


Advertising rate card available on request. Back issues and current subscriptions are also available in microform. 
The DM prices shown include postage and insurance. For subscription rates in the Americas, Japan, UK and Eire please contact your nearest 
Pergamon office. Prices and proposed publication dates are subject to change without prior notice. 

















PERGAMON PRESS 


BOOKS 
RESEARCH JOURNALS 
MAJOR REFERENCE WORKS 
CURRENT AWARENESS SERVICES 
DATABASES 
REPORTS 


PUBLISHERS 


of 
RESEARCH and LEARNING in all forms: Traditional and Electronic 


Announcing New and Newly Acquired Journals 
For 1989 & 1990 


Life Sciences & Medicine Engineering, Materials Science & Energy 
Annual Review of Addictions Research and Treatment Advanced Materials Abstracts 
Biotechnology Education Computer Optics USSR 
Cancer Communications Computing Systems in Engineering 
Cellular Signalling Mechatronics 


HEC Forum Nonwovens Abstracts 
The Journal of JASTRO World Ceramics Abstracts 
Journal of Vestibular Research Welding Abstracts 
Progress in Growth Factor Research 
Quality Assurance in Health Care 


Physical Sciences Social & Behavioural Sciences 


Astronomy Quarterly Contemporary European Affairs 
Atmospheric Environment, Part B: Urban Atmosphere European Judaism 


Quaternary International Language Sciences 
Tetrahedron: Assymetry Journal of Neurolinguistics 


Pergamon Press plc P Pergamon Press, Inc. 
Headington Hill Hall Pergamon Press Fairview Park 
Oxford OX3 OBW, UK Member of Maxwell Macmillan Pergamon Publishing Corporation Elmsford, NY 10523, USA 














JOURNAL OF THE MECHANICS AND PHYSICS OF SOLIDS 


1. GENERAL 


The Journal brings together and relates research by mathematicians, engineers, metallurgists and 
physicists on the properties of constructional materials, giving equal coverage of theory and experiment; it 
is particularly concerned with the development and application of fundamental ideas bearing on the 
connections between continuum and micro-structural properties of materials and on the solution of 
practically significant problems. 


|. Papers should be submitted to the Editor-in-Chief. Papers submitted for publication will naturally 
be judged from the standpoint of scientific originality and literary style. The preparation of progress 
reports and review articles is invited by the Editor-in-Chief from time to time and it is therefore requested 
that no uninvited contributions of the kind be submitted. The Journal appears bi-monthly, six issues 
making up one volume. 

2. Only papers that have not been previously published can be accepted and the authors must agree 
not to publish elsewhere the papers submitted to and accepted by the Journal. All requests to reproduce a 
paper either wholly or in part should be addressed to the publishers. 

3. Reprints of each paper are available to authors and a reprint order form will be sent with the 
proofs. 


Il. ScripT REQUIREMENTS 


1. Authors must submit two copies (one to be a top copy) of all manuscripts and all ancillary material. 

2. Scripts should be typed, on one side of the paper only, double spaced with good margins. The text 
must be ready for printing and should be carefully checked for errors. Authors will receive proofs for 
correction. 

3. Authors should carefully identify any hand-written symbols. It is preferred that scripts are typed 
throughout. 

4. Half-tone illustrations are to be restricted to the minimum necessary. They should accompany the 
script mounted on separate sheets. Line drawings should include all relevant details and it is particularly 
requested that originals and not photo-prints should be sent. Any lettering on diagrams—other than 
Greek letters, symbols, etc.—should be stencilled and should be sufficiently large and bold to permit 
reproduction when the diagram has been reduced to a size suitable for inclusion in the journal. 
Photographs should be enlarged sufficiently to permit their clear reproduction in half-tone. If words or 
numbers appear in photographs, two copies are requested, one clearly printed and the other without 
inscription. Tables should be of suitable quality for direct reproduction. For guidance on how best to 
prepare your tables for photographic reproduction please contact the Reprographic Section of your 
institution who will give you assistance. 

5. References to published literature should be quoted in the text as follows: SMITH (1950)—the date 
of the publication, in parentheses, following the author’s name. References should be listed together at the 
end of each paper and not given as footnotes. The title of the journal should be abbreviated as in World 
List of Scientific Periodicals (4th edn). They should be arranged in alphabetical order (first author’s sur- 
name) to appear as follows: 

KENNEDY,A. J. 1953 J. Mech. Phys. Solids 1, 172. 
Roserts, G. H. 1952 Chem. Engng Sci.1, 196. 


6. The text should be concise and in a readily understandable style. The technical description of the 
methods used should only be given in detail when such methods are new. The essential contents of each 
paper should be briefly recapitulated in a summary on the title page. 

7. To conserve space authors are requested to mark less important portions of the paper (such as 
descriptions of methods, record of experimental results, etc.) for printing in smaller type. 

8. Corrections to the proofs should be restricted to the printer’s errors only. Any other substantial 
changes at this stage may be charged to the author. 

9. It would be appreciated if authors would notify the Publisher of any change of address which 
occurs whilst their papers are in the process of publication. 

10. The manuscript and diagrams will be discarded one month after publication unless the Publisher 
is requested to return original material to the author. 








JOURNAL OF THE MECHANICS AND PHYSICS OF 
SOLIDS 


VOLUME 38, NUMBER 1 1990 


CONTENTS 


I. F. COLLINS 


N. P. Kruyt 

J.O. CHUNG, Jin YU and 
S. H. HONG 

R. L. WEAVER 

W. W. GERBERICH, 

D. L. DAvipsoNn and 


M. KACZOROWSKI 


W.KREHER 


Corrigendum 


Announcement 


Indexed/Abstracted in: 


Plane strain characteristics theory for soils and 
granular materials with density dependent yield 
criteria 


An analysis of the generalized double-sliding 
models for cohesionless granular materials 


Steady-state creep crack growth by continually 
nucleating cavities 


Diffusivity of ultrasound in polycrystals 
Experimental and theoretical strain distributions 


for stationary and growing cracks 


Residual stresses and stored elastic energy of 
composites and polycrystals 


International Centre for Mechanical Sciences, 
Preliminary Programme 1990 


Current Contents, INSPEC, Appl. Mech Rey., Cambridge Scient. Abstr., ISSN 0022-5096 
Mechanics, Math. Rev. and Pascal/CNRS Database JMPSAS8 38 (1) 1-132 (1990) 


Published by 


PERGAMON PRESS OXFORD * NEWYORK 
BEIJING - FRANKFURT - SAO PAULO - SYDNEY - TOKYO - TORONTO 


Printed in Great Britain at Aberdeen University Press 220 








Reproduced with the permission of Pergamon Press Inc., by University 
Microfilms Inc. Duplication or resale without permission is prohibited. 





